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Post-Newtonian celestial dynamics is a relativistic theory of motion of massive bodies and test 
particles under the influence of relatively weak gravitational forces. Standard approach for develop- 
ment of this theory relies upon the key concept of the isolated astronomical system supplemented 
by the assumption that the background space-time is flat. The standard post-Newtonian theory 
of motion was instrumental in explanation of the existing experimental data on binary pulsars, 
satellite and lunar laser ranging, and in building precise ephemerides of planets in the solar system. 
Recent studies of the formation of large-scale structure in our universe indicate that the standard 
post-Newtonian mechanics fails to describe more subtle dynamical effects in motion of the bodies 
comprising the astronomical systems of larger size - galaxies and clusters of galaxies - where the 
Riemann curvature of the expanding FLRW universe interacts with the local gravitational field of 
the astronomical system and, as such, can not be ignored. 

The present paper outlines theoretical principles of the post-Newtonian mechanics in the expand- 
ing universe. It is based upon the gauge-invariant theory of the Lagrangian perturbations of cosmo- 
logical manifold caused by an isolated astronomical N-body system (the solar system, a binary star, 
a galaxy, a cluster of galaxies) . We postulate that the geometric properties of the background man- 
ifold are described by a homogeneous and isotropic Friedman- Lemaitre-Robertson- Walker (FLRW) 
metric governed by two primary components - the dark matter and the dark energy. The dark 
matter is treated as an ideal fluid with the Lagrangian taken in the form of pressure along with the 
scalar Clebsch potential as a dynamic variable. The dark energy is associated with a single scalar 
field with a potential which is hold unspecified as long as the theory permits. Both the Lagrangians 
of the dark matter and the scalar field are formulated in terms of the field variables which play a 
role of generalized coordinates in the Lagrangian formalism. It allows us to implement the powerful 
methods of variational calculus to derive the gauge-invariant field equations of the post-Newtonian 
celestial mechanics of an isolated astronomical system in an expanding universe. These equations 
generalize the field equations of the post-Newtonian theory in asymptotically-flat spacetime by tak- 
ing into account the cosmological effects explicitly and in a self-consistent manner without assuming 
the principle of liner superposition of the fields or a vacuole model of the isolated system, etc. The 
field equations for matter dynamic variables and gravitational field perturbations are coupled in the 
most general case of arbitrary equation of state of matter of the background universe. 

We introduce a new cosmological gauge which generalizes the de Donder (harmonic) gauge of 
the post-Newtonian theory in asymptotically flat spacetime. This gauge significantly simplifies the 
gravitational field equations and allows to find out the approximations where the field equations 
can be fully decoupled and solved analytically. The residual gauge freedom is explored and the 
residula gauge transformations are formulated in the form of the wave equations for the gauge 
functions. We demonstrate how the cosmological effects interfere with the local system and affect 
the local distribution of matter of the isolated system and its orbital dynamics. Finally, we worked 
out the precise mathematical definition of the Newtonian limit for an isolated system residing on 
the cosmological manifold. The results of the present paper can be useful in the solar system for 
calculating more precise ephemerides of the solar system bodies on extremely long time intervals, in 
galactic astronomy to study the dynamics of clusters of galaxies, and in gravitational wave astronomy 
for discussing the impact of cosmology on generation and propagation of gravitational waves emitted 
by coalescing binaries and/or merging galactic nuclei. 
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I. NOTATIONS 

This section summarizes notations used in the present paper. We use G to denote the universal gravitational 
constant and c for the ultimate speed in Minkowski spacetime. Every time, when there is no confusion about the 
system of units, we shall choose a geometrized system of units such that G = c = 1. We put a bar over any function 
that belongs to the background manifold of the FLRW cosmological model. Any function without such a bar belongs 
to the perturbed manifold. 

The notations used in the present paper are as follows: 

• Greek indices a, f3, 7, . . . run through values 0, 1, 2, 3, and Roman indices k, . . . take values 1, 2, 3, 

• Einstein summation rule is applied for repeated (dummy) indices, for example, P a Q a = P°Q + P 1 Qi+ P 2 Q 2 + 
P 3 Q 3 , and P'Q l = P 1 Q 1 +P 2 Q 2 + P 3 Q 3} 

• g a f) is a full metric on the cosmological spacetime manifold, 

• g a /3 is the FLRW metric on the background spacetime manifold, 

• f a /3 is the metric on the conformal spacetime manifold, 

• Va/3 = diag{— 1, +1, +1, +1} is the Minkowski metric, 

• T and X % = {A, Y, Z} are the coordinate time and isotropic spatial coordinates on the background manifold, 

• X a = {A , A 4 } = {c?7, A 4 } are the conformal coordinates with 77 being a conformal time, 

• x a — {x° 7 x 1 } = {ct,x 1 } is an arbitrary coordinate chart on the background manifold, 

• a bar, F above a geometric object F, denotes the unperturbed value of F on the background manifold, 

• a prime F' = dF/dr/ denotes a total derivative with respect to the conformal time 77, 

• a dot F = dF/dij denotes a total derivative with respect to the cosmic time T, 

• d a = d/dx a is a partial derivative with respect to the coordinate x a 7 

• a comma with a following index F Q = d a F is another designation of a partial derivative with respect to a 
coordinate x a , 

• a vertical bar, Fi a denotes a covariant derivative of a geometric object F (a scalar, a vector, a tensor) with 
respect to the background metric g a /3, 

• a semicolon, F. a denotes a covariant derivative of a geometric object F (a scalar, a vector, a tensor) with respect 
to the conformal metric f a /3, 

• the tensor indices of geometric objects on the background manifold are raised and lowered with the background 
metric g a /3, 

• the tensor indices of geometric objects on the conformal spacetime are raised and lowered with the conformal 
metric f Q/ 3, 

• the scale factor of the FLRW metric is denoted as R = R(T), or as a = a(rf) = R[T(rj)}, 

• the Hubble parameter, H= R/R 7 and the conformal Hubble parameter, 'K = a' /a. 
Other notations will be introduced and explained in the main text of the paper. 
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II. INTRODUCTION 



Post-Newtonian celestial mechanics is a branch of fundamental gravitational physics [12, 58, 92] that deals with the 
theoretical concepts and experimental methods of measuring gravitational fields and testing general theory of relativity 
both in the solar system and beyond [97, 104]. In particular, the relativistic celestial mechanics of binary pulsars 
(see [72], and references therein) was instrumental in providing conclusive evidence for the existence of gravitational 
radiation as predicted by Einstein's theory [93, 101]. 

Over the last few decades, various groups within the International Astronomical Union (IAU) have been active 
in exploring the application of the general theory of relativity to the modeling and interpretation of high-accuracy 
astrometric observations in the solar system and beyond. A Working Group on Relativity in Celestial Mechanics and 
Astrometry was formed in 1994 to define and implement a relativistic theory of reference frames and time scales. This 
task was successfully completed with the adoption of a series of resolutions on astronomical reference systems, time 
scales, and Earth rotation models by 24th General Assembly of the IAU, held in Manchester, UK, in 2000. The IAU 
resolutions are based on the first post-Newtonian approximation of general relativity which is a conceptual basis of 
the fundamental astronomy in the solar system [91]. 

The mathematical formalism of the post-Newtonian approximations is getting progressively complicated as one 
goes from the Newtonian to higher orders [23, 88]. For this reason the theory has been primarily developed under 
a basic assumption that the background spacetime is asymptotically flat. Mathematically, it means that the full 
spacetime metric, g a p, is decomposed around the background Minkowskian metric, r\ a $ — diag(— 1, 1, 1, 1), into a 
linear combination 



9af) = Val3 + h a p , (1) 

where the perturbation, 

h a ^c-*h%+c-^l + c-^l + ... , (2) 

is the post- Newtonian series with respect to the powers of 1/c, where c is the ultimate speed in general relativity 
(equal to the speed of light). Post-Newtonian approximations is the method to determine h a p by solving Einstein's 
field equations with the tensor of energy-momentum of matter of a localized astronomical system, T Q( g($, h a p), taken 
as a source of the gravitational field h a p, by iterations starting from h a p = in the expression for % a p. The 
solution of the field equations and the equations of motion of the astronomical bodies are derived in some coordinates 
r a = {ct, r} where t is the coordinate time, and r = {x, y, z} are spatial coordinates. The post-Newtonian theory in 
asymptotically-flat spacetime has a well-defined Newtonian limit determined by: 

1) equation for the Newtonian potential, U — 



Jv \i - i I 

where p = c~ 2 1 00 , is the density of matter producing the gravitational field, 

2) equation of motion for massive particles 

r = VU , (4) 

where V = {d x ,d y ,d z } is the operator of gradient, r = r(t) is time-dependent position of a particle (worldline 
of the particle), and the dot denotes a total derivative with respect to time t, 

3) equations of motion for light (massless particles) 

f = . (5) 



These equations are considered as fundamentals for creation of astronomical ephemerides of celestial bodies in the 
solar system [58] and in any other localized system of self-gravitating bodies like a binary pulsar [72]. In all practical 
cases they have to be extended to take into account the post-Newtonian corrections sometimes up to the 3-d post- 
Newtonian order of magnitude [105]. It is important to notice that in the Newtonian limit the coordinate time t of the 
gravitational equations of motion (4), (5) coincides with the proper time of observer r that is practically measured with 
an atomic clock. The formalism of the present paper has been employed in [60] to check the theoretical consistency of 
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(3)-(5) and to analyse the outcome of some experiments like the anomalous Doppler effect discovered by J. Anderson 
et al [3, 4] in the orbital motion of Pioneer 10 and 11 space probes. 

So far, the post-Newtonian theory was mathematically successful ("unreasonably effective" as Clifford Will states 
[105]) and passed through numerous experimental tests with a flying color. Nevertheless, it hides several pitfalls. 
The first one is the problem of convergence of the post-Newtonian series and regularization of divergent integrals 
that appear in the post- Newtonian calculations at higher post- Newtonian orders [88]. The second problem is that 
the background manifold is not asymptotically- flat Minkowskian spacetime but the FLRW metric, g a p. We live in 
the expanding universe which rate of expansion is determined by the Hubble constant H n . Therefore, the right thing 
would be to replace the post- Newtonian decomposition (1) with a more adequate post-Friedmannian series [96] 

g a /3 = 9af3 + X a [ 3 , (6) 

where 

= *$ } + + H 2 k% } + ..., (7) 

is the metric perturbation around the cosmological background represented as a series with respect to the Hubble 
parameter, H. Each term of the series has its own expansion into post- Newtonian series like (2). Generalization 
of the theory of post-Newtonian approximations from the Minkowski spacetime to that of the expanding universe is 
important for extending the applicability of the post-Newtonian celestial mechanics to testing cosmological effects and 
for more deep understanding of the process of formation of the large-scale structure in the universe and gravitational 
interaction between galaxies and clusters of galaxies. 

Whether cosmological expansion affects gravitational dynamics of bodies inside a localized astronomical system was 
a matter of considerable efforts of many researchers [11, 30, 31, 62, 77, 78, 89]. A recent article [16] summarizes the 
previous results and provides the reader with a number of other valuable resources. Most of the previous works on 
celestial mechanics in cosmology were based on assumption of spherical symmetry of gravitational field and matching 
two (for example, Schwarzschild and Friedmann) exact solutions of Einstein's equations. The matching was achieved 
in many different ways. McVittie's solution [78] is perhaps the most successful mathematically but yet lacks a clear 
physical interpretation [16]. Moreover, its practical application is doubtful since it is valid only for spherically- 
symmetric case. 

We need a precise mathematical formulation of the post-Newtonian theory for a self-gravitating localized astro- 
nomical system not limited by the assumption of the spherical symmetry, embedded to the expanding universe and 
coupled through the gravitational interaction with the time-dependent background geometry. Theoretical description 
of the post-Newtonian theory for a localized astronomical system in expanding universe should correspond in the limit 
of vanishing H to the post-Newtonian theory obtained in the asymptotically-flat spacetime. Such a description will 
allow us to directly compare the equations of the standard post-Newtonian celestial mechanics with its cosmological 
counterpart. Therefore, the task is to derive a set of the post-Newtonian equations in cosmology in some coordi- 
nates introduced on the background manifold, and to map them onto the set of the Newtonian equations (3)-(5) 
in asymptotically-flat spacetime. Such a theory of the post-Newtonian celestial mechanics would be of a paramount 
importance for extending the tools of experimental gravitational physics to the field of cosmology, for example, to 
properly formulate the cosmological extension of the PPN formalism [103]. The present article discusses the main ideas 
and principal results of such a theoretical approach in the linearized approximation with respect to the gravitational 
perturbations of the cosmological background caused by the presence of a localized astronomical system. 

The paper is organized as follows. In the following section we describe a brief history of the development of the 
theory of cosmological perturbations. Section IV describes the Lagrangian of gravitational field, the matter of the 
background cosmological model, and an isolated astronomical system which perturbs the background cosmological 
manifold. Section V describes the geometric structure of the background spacetime manifold of the cosmological 
model and the corresponding equations of motion of the matter and field variables. Section VI introduces the reader 
to the theory of the Lagrangian perturbations of the cosmological manifold and the dynamic variables. Section VII 
makes use of the preceding sections in order to derive the field equations in the gauge-invariant form. Beginning 
from section VIII we focus on the spatially-flat universe in order to derive the post-Newtonian field equations that 
generalize the post-Newtonian equations in the asymptotically-flat spacetime. These equations are coupled in the 
scalar sector of the proposed theory. Therefore, we consider in section IX a few particular cases when the equations 
can be fully decoupled one from another, and solved in terms of the retarded potentials. Appendix provides a proof of 
the Lorentz-invariance of the retarded potentials for the wave equations describing propagation of weak gravitational 
and sound waves on the background cosmological manifold. 
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III. BRIEF HISTORY OF COSMOLOGICAL PERTURBATION THEORY 

In order to solve the problem of the interaction of the gravitational field of an isolated astronomical system with 
expanding universe one has to resort to the theory of cosmological perturbations. The immediate goal of cosmological 
perturbation theory is to relate the physics of the early universe to CMB anisotropy and large-scale structure and to 
provide the initial conditions for numerical simulations of structure formations. The ultimate goal of this theory is 
to establish a mathematical link between the fundamental physical laws at the Planck epoch and the output of the 
gravitational wave detectors which are the only experimental devices being able to map parameters of the universe at 
that time [56]. 

Originally, two basic approximation schemes for calculation of cosmological perturbations have been proposed by 
Lifshitz with his collaborators [68, 69] and, later on, by Bardeen [7]. Lifshitz [68] worked out a coordinate-dependent 
theory of cosmological perturbations while Bardeen [7] concentrated on finding the gauge-invariant combinations for 
perturbed quantities and derivation of a perturbation technique based on gauge-invariant field equations. At the 
same time, Lukash [73] had suggested an original approach for deriving the gauge-invariant scalar equations based on 
the thermodynamic theory of the Clebsch potential also known in cosmology as the scalar velocity potential [64, 90] 
or the Taub potential [90, 95]. It turns out that the variational principle with a Lagrangian of cosmological matter 
formulated in terms of the Clebsch potential, is the most useful mathematical device for developing the theory of 
relativistic celestial mechanics of localized astronomical systems embedded in expanding cosmological manifold. It is 
for this reason, we use the Clebsch potential in the present paper. 

A few words of clarification regarding a comparison between Lifshitz's and Bardeen's approaches should be relayed 
to the reader. The approach established in the papers by Lifshitz [68, 69], is fully correct. Lifshitz decided to work 
in synchronous gauge and realized that this fixing of coordinates allows for a residual gauge ambiguity, which can 
also be fixed by picking a synchronous, comoving coordinate system. Bardeen [7] wrote his paper, because there was 
some confusion in the 1970s about that issue (which could have beed avoided if people would have studied Lifshitz's 
papers carefully). He demonstrated, that any coordinate could be chosen and that there exist quantities which are 
independent of that choice, which he identified with the physical degrees of freedom. However, this is not where the 
story of the cosmological perturbation theory ends. Closer inspection shows that what is really relevant is not the 
choice of coordinates (which do not have a physical meaning), rather the choice of spacetime foliation is relevant, for 
example, it makes a physical difference if one defines the Harrison-Zel'dovich spectrum [46, 107] of the primordial 
perturbations on a synchronous, comoving hypersurface or a shear free hypersurface. Pitfalls in understanding this 
issue are subtle and, sometimes, may be not easily recognized (see [25, 42, 76] and [40] for a detailed discussion of the 
role of foliations in cosmology and in general relativity) . 

In the years that followed, the gauge-invariant formalism was refined and improved by Durrer and Straumann 
[27, 28], Ellis et al. [32-34] and, especially, by Mukhanov et al. [81, 82]. Irrespectively of the approach a specific 
gauge must be fixed in order to solve equations for cosmological perturbations. Any gauge is allowed and its particular 
choice is simply a matter of convenience. Imposing a gauge condition eliminates four gauge degrees of freedom in the 
cosmological pertrubations and brings the differential equations for them to a solvable form. Nonetheless, the residual 
gauge freedom originated from the tensor nature of the gravitational field remains. This residual gauge freedom leads 
to appearance of unphysical perturbations which must be disentangled from the physical modes. Lifshitz theory of 
cosmological perturbations [68, 69] is worked out in a synchronous gauge and contains the spurious modes but they 
are easily isolated from the physical perturbations [41]. Other gauges used in cosmology are described in Bardeen's 
paper [7] and used in cosmological perturbation theory. Among them, the longitudinal (conformal or Newtonian) 
gauge is one of the most common. This gauge is advocated by Mukhanov [81] because it removes spurious coordinate 
degrees of freedom out of scalar perturbations. Detailed comparison of the cosmological perturbation theory in the 
synchronous and conformal gauges was given by Ma and Bertschinger [74]. 

Unfortunately, none of the previously known cosmological gauges can be applied for analysis of the cosmological 
perturbations caused by localized matter distributions like an isolated astronomical system which can be a single 
star, a planetary system, a galaxy, or even a cluster of galaxies. The reason is that the synchronous gauge has no 
the Newtonian limit and is applicable only for freely falling test particles while the longitudinal gauge separates the 
scalar, vector and tensor modes present in the metric tensor perturbation in the way that is incompatible with the 
technique of the approximation schemes having been worked out in asymptotically flat space-time [58] . We also notice 
that standard cosmological perturbation technique often operates with harmonic (Fourier) decomposition of both the 
metric tensor and matter perturbations when one is interested in statistical statements based on the cosmological 
principle. This technique is unsuitable and must be avoided in sub-horizon approximation for working out the 
post-Newtonian celestial mechanics of self-gravitating isolated systems. Current paradigm is that the cosmological 
generalization of the Newtonian field equations of an isolated gravitating system like the solar system or a galaxy 
or a cluster of galaxies can be easily obtained by just making use of the linear principle of superposition with a 
simple algebraic addition of the local system to the tensor of energy momentum of the background matter. It is 
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assumed that the superposition procedure is equivalent to operating with the Newtonian equations of motion derived 
in asymptotically-flat spacetime and adding to them ("by hands") the tidal force due to the presence of the external 
universe (see, for example, [77]). Though such a procedure may look pretty obvious it lacks a rigorous mathematical 
analysis of the perturbations induced on the background cosmological manifold by the local system. This analysis 
should be done in the way that embeds cosmological variables to the field equations of standard post-Newtonian 
approximations not by "hands" but by precise mathematical technique which is the goal of the present paper. The 
variational calculus on manifolds is the most convenient for joining the standard theory of cosmological perturbations 
with the post-Newtonian approximations in asymptotically-flat spacetime. It allows us to track down the rich interplay 
between the perturbations of the background manifold with the dynamic variables of the local system which cause 
these perturbations. The output is the system of the post-Newtonian field equations with the cosmological effects 
incorporated to them in a physically-transparent and mathematically-rigorous way. This system can be used to solve 
a variety of physical problems starting from celestial mechanics of localised systems in cosmology to gravitational wave 
astronomy in expanding universe that can be useful for deeper exploration on scientific capability of such missions as 
LISA and Big Bang Observer (BBO) [22] 

In fact, the problem of whether the cosmological expansion affects the long-term evolution of an isolated N-body 
system (galaxy, solar system, binary system, etc.) has a long controversial history. The reason is that there was 
no an adequate mathematical formalism for describing cosmological perturbations caused by the isolated system so 
that different authors have arrived to opposite opinions. It seems that McVittie [78] was first who had considered 
the influence of the expansion of the universe on the dynamics of test particles orbiting around a massive point-like 
body immersed to the cosmological background. He found an exact solution of the Einstein equations in his model 
which assumed that the mass of the central body is not constant but decreases as the universe expands. Einstein and 
Straus [30, 31] suggested a different approach to discuss motion of particles in gravitationally self-interacting systems 
residing on the expanding background. They showed that a Schwarzschild solution could be smoothly matched to the 
Friedman universe on a spherical surface separating the two solutions. Inside the surace ("vacuole") the motion of 
the test particles is totally unaffected by the expansion. Thus, Einstein and Straus [30, 31] concluded that the cosmic 
expansion is irrelevant for the Solar system. Bonnor [11] generalized the Einstein-Straus vacuole and matched the 
Schwarzschild region to the inhomogeneous Lemaitre-Tolman-Bondi model thus, making the average energy density 
inside the vacuole be independent of the exterior energy density while in the Einstein-Straus model they must be 
equal. Bonnor [11] concluded that the local systems expand but at a rate which is negligible compared with the 
general cosmic expansion. Similar conclusion was reached by Mashhoon et al. [77] who analysed the tidal dynamics 
of test particles in the Fermi coordinates. 

The vacuole solutions are not appropriate for adequate physical solution of the N-body problem in the expanding 
universe. There are several reasons for it. First, the vacuole is spherically-symmetric while majority of real astro- 
nomical systems are not. Second, the vacuole solution imposes physically unrealistic boundary conditions on the 
matching surface that relates the central mass to the size of the surface and to the cosmic energy density. Third, the 
vacuole is unstable against small perturbations. In order to overcome these difficulties a realistic approach based on 
the approximate analytic solution of the Einstein equations for the N-body problem immersed to the cosmological 
background, is required. In the case of a flat space-time there are two the most advanced techniques for finding 
approximate solution of the Einstein equations describing gravitational field of an isolated astronomical system. The 
first is called the post-Newtonian approximations and we have briefly discussed this technique in the introduction. 
The second technique is called post-Minkowskian approximations [23]. The post-Newtonian approximations is appli- 
cable to the systems with weak gravitational field and slow motion of matter. The post-Minkowskian approximations 
also assume that the field is weak but does not imply any limitation on the speed of matter. The post-Newtonian 
iterations are based on solving the elliptic-type Poisson equations while the post-Minkowskain approach operates with 
the hyperbolic-type Dalambert equations. The post-Minkowskian approximations naturally includes description of 
the gravitational radiation emitted by the isolated system while the post-Newtonian scheme has to use additional 
mathematical methods to describe generation of the gravitational waves [17]. In the present paper we concentrate 
on the development of a generic scheme for calculation of cosmological perturbations caused by a localized distribu- 
tion of matter which preserves many advantages of the post-Minkowskian approximation scheme. The cosmological 
post-Newtonian approximations are derived from the general perturbation scheme by making use of the slow-motion 
expansion with respect to a small parameter v/c where v is the characteristic velocity of matter in the N-body system 
and c is the fundamental speed. 

There were several attempts to work out a physically-adequate and mathematically-rigorous approximation schemes 
in general relativity in order to construct and to adequately describe small-scale self-gravitating structures in the 
universe. The most notable works in this direction have been done by Kurskov and Ozernoi [63], Futamase et al. 
[10, 37, 38], Buchert and Ehlers [14, 29], Mukhanov et al. [1, 81-83], Zalaletdinov [106]. These approximation schemes 
have been designed to track the temporal evolution of the cosmological perturbations from a very large down to a small 
scale up to the epoch when the perturbation becomes isolated from the expanding cosmological background. These 
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approaches looked hardly connected between each other until recent works by Clarkson et al [20, 21], Li & Schwarz 
[66, 67], Rasanen [87], Buchert & Rasanen [15] and Wiegand & Schwarz [102]. In particular, Wiegand & Schwarz 
[102] have shown that the idea of cosmic variance (that isa standard way of thinking) is closely related to the cosmic 
averages defined by Buchert and Ehlers [14, 29]. All researchers agree that the second and higher-order non-linear 
approximations are important to understand the back-reaction of the cosmological perturbations propagating on the 
cosmological background used in the linearized theory (see, for example, [1, 37, 50, 51, 82, 106]). 

Papers [24, 61, 86] attempted to construct an approximation scheme being compatible with both the post-Newtonian 
and post-Minkowskian approximations in asymptotically-flat space-time and the gauge-invariant theory of cosmolog- 
ical perturbations caused by a localized astronomical system. We have succeeded in solving this problem in the work 
[61] which assumes the dust-dominated background universe with spatial curvature k — 0. We remind that in standard 
Bardeen's approach [7] the metric tensor perturbations h a p are decomposed in irreducible scalar, vector, and tensor 
parts which are combined with themselves and with matter perturbations in order to obtain some gauge-invariant 
quantities that do not contain spurious modes invoked by the freedom in doing coordinate transformations on cos- 
mological background manifold. Bardeen [7] reformulates Einstein's field equations in terms of these gauge-invariant 
variables which are further decomposed in Fourier harmonics. The field equations become then of the Helmholtz 
type and are solved by constructing Green's function. This specific procedure of Bardeen's approach is incompatible 
with the post-Newtonian or post-Minkowskian approximations which do not decompose the metric tensor in scalar, 
vector, and tensor harmonics and do not expand them to the Fourier series. Therefore, we have used a different 
procedure based on introduction of auxiliary scalar, <f>, and vector, £, fields which are used along with the metric 
tensor perturbation h a p as basic elements for decomposing the perturbed stress-energy tensor of matter 6T a p and 
selecting from this decomposition that part of 5T a p which has the same transformation property as the perturbed 
Einstein tensor SG a /3. This process makes the rest of the perturbation of 8T a p gauge- invariant so that it can be 
identified with the bare (external) perturbation imposed on the cosmological background by the presence of a material 
system like a single star, solar system, galaxy, etc. The auxiliary scalar, (p, and vector, £, fields are also determined 
by this procedure. For example, it is proved out [61] that the vector field £ is identically equal to zero while the 
scalar field (f> is found from the equation following from the Bianchi identity of the perturbed Einstein equations. The 
entire approach is gauge-invariant but the equations for the scalar and metric perturbations are strongly coupled in 
general case. We have shown that there is a special cosmological gauge generalizing the harmonic gauge of general 
relativity in such a way that the reduced field equations are completely decoupled and significantly simplified. More 
specifically, the linearized Einstein equations for the metric tensor perturbations hoi and hij are decoupled both from 
each other and from hoo component which couples only with the auxiliary scalar field <f>. However, it turns out that 
our gauge [61] admits a simple linear combination, x, of hoo, hkk, an d <ft such that the equation for \ decouples from 
any other perturbation. The equations for x, hoi and hij are of a wave-type and have the bare stress energy-tensor 
of matter as a source in their right-hand sides. These equations have simple Green functions given in terms of the 
retarded integrals with the Minkowskian null cone defined by the conformally-flat part of the FLRW metric. 

In the work [86] this linearized approach has been extended to the background cosmological models governed by 
a perfect fluid with the barotropic equation of state p — qe, where p and e are pressure and energy density of the 
background matter respectively, and q is a constant parameter taking value in the range from -1 to +1. We have 
shown [86] that the overall perturbative scheme for calculation of the cosmological perturbations in such model can 
be significantly streamlined and simplified if one formally replaces the stress-energy tensor of the perfect fluid with 
one of a classic scalar field minimally coupled with metric. A specific (exponential) form of the potential V($) of the 
scalar field $ is fixed by two conditions: 

1. it reproduces all functional relationships of the background FLRW cosmological model; 

2. it maintains the background equation of state p = qe. 

Although a minimally coupled scalar field can be viewed to a certain extent as a perfect fluid one should keep in 
mind that its barotropic equation of state does not hold generally in the perturbed universe (see [13] and discussion 
in [42, 76]). This may impose some technical difficulties in handling the mathematical analysis of cosmological 
perturbations caused by localized distributions of matter. We explain how to get around these difficulties by making 
use of the Lagrangian-based variational technique of the gauge-invariant perturbations of curved manifolds. This 
makes our perturbative approach [61, 86] more efficient in developing the post-Newtonian celestial mechanics in 
cosmology as compared with the standard technique [7, 27, 28, 32-34, 68, 69, 74, 81, 82]. 

Development of observational cosmology and gravitational wave astronomy demands to extend the linearized theory 
of cosmological perturbations to second and higher orders of approximation. A fair number of works have been devoted 
to solving this problem. Non-linear perturbations of the metric tensor and matter affect evolution of the universe and 
this back-reaction of the perturbations should be taken into account. This requires derivation of the effective stress- 
energy tensor for cosmological perturbations like freely-propagating gravitational waves and scalar field [1, 81-83]. 
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The laws of conservation for the effective stress-energy tensor are important for better understanding of physics of 
the expanding universe [6, 43]. 

In the present paper we construct a non-linear theory of cosmological perturbations for isolated systems which 
generalizes the post-Minkowskian approximation scheme for calculation perturbations of gravitational field in asymp- 
totically flat space-time. We rely upon the basic results of the linearized gauge-invariant theory from our previous 
works [61, 86] in order to derive a decoupled system of equations for quadratic cosmological perturbations of a spatially 
flat FLRW universe. We implement the Lagrangian-based theory of dynamical perturbations of gravitational field 
on a curved background manifold which has been worked out in [44, 85] (see also [6]). This theory has a number of 
specific advantages over other perturbation methods among which the most important are: 

• Lagrangian-based approach is covariant and can be implemented for any curved background spacetime that is 
an exact solution of the Einstein gravity field equations; 

• the system of the partial differential equations describing dynamics of the perturbations is determined by a 
dynamic Lagrangian L D which is derived from the total Lagrangian L by making use of its Taylor expansion 
with respect to the perturbations and accounting for the background field equations. The dynamic Lagrangian 
L D defines the conserved quantities for the perturbations (energy, angular momentum, etc.) that depend on 
the symmetries of the background manifold; 

• the dynamic Lagrangian L D and the corresponding field equations for the perturbations are gauge- invariant in 
any order of the perturbation theory. Gauge transformations map the background manifold onto itself and are 
associated with arbitrary (analytic) coordinate transformations on the background space-time; 

• the entire perturbation theory is self-reproductive and is extended to the next perturbative order out of a 
previous iteration so that the linearized approximation is the basic starting point. 



We accept the Einstein's theory of general relativity and consider a universe filled up with matter consisting of three 
components. The first two components are: (1) an ideal fluid composed of particles of one type with transmutations 
excluded; (2) a scalar field; and (3) a matter of the localized astronomical system. The ideal fluid consists of baryons 
and cold dark matter, while the scalar field describes dark energy [2]. We assume that these two components do not 
interact with each other directly, and are the source of the Friedmann-Lemitre-Robertson- Walker (FLRW) geometry. 
There is no dissipation in the ideal fluid and in the scalar field so that they can only interact through the gravitational 
field. It means that the equations of motion for the fluid and the scalar field are decoupled, and we can calculate 
their evolution separately. In other words, the Lagrangian of the ideal fluid and that of the scalar field depend only 
on their own field variables. 

The tensor of energy-momentum of matter of the localized astronomical system is not specified in agreement with 
the approach adopted in the post-Newtonian approximation scheme developed in the asymptotically-flat spacetime 
[23, 57]. This allows us to generate all possible types of cosmological perturbations: scalar, vector and tensor modes. 
We are the most interested in developing our formalism for application to the astronomical system of massive bodies 
bound together by intrinsic gravitational forces like the solar system, galaxy, or a cluster of galaxies. It means that 
our approach admits a large density contrast between the background matter and the matter of the localized system. 
The localized system perturbs the background matter and gravitational field of FLRW universe locally but it is 
not included to the matter source of the background geometry, at least, in the approximation being linearized with 
respect to the metric tensor perturbation. Our goal is to study how the perturbations of the background matter and 
gravitational field are incorporated to the gravitational field perturbations of the standard post-Newtonian theory of 
relativistic celestial mechanics. 

Let us now consider the action functional and the Lagrangian of each component. 



IV. LAGRANGIAN AND FIELD VARIABLES 



A. 



The Action Functional 



We shall consider a theory with the action functional 




(8) 
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where the integration is performed over the entire spacetime manifold M. The Lagrangian £ is comprised of four 
terms 

L = £ s + £ m + £ q + £ p , (9) 

where £ s , £ m , £ q are the Lagrangians of gravitational field, the dark matter, the scalar field that governs the acceler- 
ated expansion of the universe [39], and L p is the Lagrangian describing the source of the cosmological perturbations. 
Gravitational field Lagrangian is 

£ s = --±-V=9R, (10) 

WTT 

where R is the Ricci scalar built of the metric g a p and its first and second derivatives [79]. Other Lagrangians depend 
on the metric and the matter variables. Correct choice of the matter variables is a key element in the development 
of the Lagrangian theory of the post-Newtonian perturbations of the cosmological manifold caused by a localized 
astronomical system. 



B. The Lagrangian of the Ideal Fluid 



The ideal fluid is characterized by the following thermodynamic parameters: the rest- mass density p m , the specific 
internal energy II m (per unit of mass), pressure p m , and entropy s m where the sub- index 'm' stands for 'matter'. 
We shall assume that the entropy of the ideal fluid remains constant, that excludes it from further consideration. 
The standard approach to the theory of cosmological perturbations preassumes that the constant entropy excludes 
rotational (vector) perturbations of the fluid component from the start, and only scalar (adiabatic) perturbations are 
generated [2, 81, 99, 100]. However, the present paper deals with the cosmological perturbations that are generated 
by a localized astronomical system which is described by its own Lagrangian (see section IV D) which is left as general 
as possible. This leads to the tensor of energy-momentum of the matter of the localized system that incorporates 
the rotational motion of matter which is the source of the rotational perturbations of the background ideal fluid. 
This extrapolates the concept of the gravitomagnetic field of the post- Newtonian dynamics of localized systems in the 
asymptotically-flat spacetime [12, 19, 58] to cosmology. Further details regarding the vector perturbations are given 
in section VII E of the present paper. 

The total energy density of the fluid 

e m -Pm(l + n m ). (11) 
One more thermodynamic parameter is the specific enthalpy of the fluid defined as 

Mm = = 1 + n m H . (12) 

Pm Pm 

In the most general case, the thermodynamic equation of state of the fluid is given by equation p m = p m (pm,n m ), 
where the specific internal energy IT m is related to pressure by the first law of thermodynamics. 
Since the entropy has been assumed to be constant, the first law of thermodynamics reads 

dn m +ft n d^^ =0. (13) 

It can be used to derive the following thermodynamic relationships 

dpm = Pmdp.m , (14) 
de m = Pmdpm , (15) 

which means that all thermodynamic quantities are solely functions of the specific enthalpy /z m , for example, p m = 
Pm(Mm), LT m = n m (/i m ), etc. The equation of state is also a function of the variable /U m , that is 

Pm = Pm(Mm) • (16) 

Derivatives of the thermodynamic quantities with respect to /j, m can be calculated by making use of equations (14) 
and (15), and the definition of the (adiabatic) speed of sound c s of the fluid 

^ = 4 > 

<9e m c 2 



10 



where the partial derivative is taken under a condition that the entropy, s m , of the fluid does not change. Then, the 
derivatives of the thermodynamic quantities take on the following form 

dp m de m c 2 dp 

where all partial derivatives are performed under the same condition of constant entropy. 

The Lagrangian of the ideal fluid is usually taken in the form of the total energy density, L m — \f zr ge m [79]. 
However, this form is less convenient for applying the variational calculus on manifolds. The above thermodynamic 
relationships and the integration by parts of the action (8) allows us to recast the Lagrangian L m = y^jtm to the 
form of pressure, L m = —\J~ r gp mi so that the Lagrangian density becomes (see [58, pp. 334-335 ] for more detail) 

£ m = -V-gPm = V^ffOm - PmPm) ■ (19) 

Theoretical description of the ideal fluid as a dynamic system on space-time manifold is given the most conveniently 
in terms of the Clebsch potential, <I> which is also called the velocity potential [90]. In the case of a single-component 
ideal fluid the Clebsch potential is introduced by the following relationship 

Mm^a = a ■ (20) 

In fact, equation (20) is a solution of relativistic equations of motion of the ideal fluid [64]. 

The Clebsch potential is a primary field variable in the Lagrangian description of the isentropic ideal fluid. The 
four-velocity is normalized to w a w a — g a pw a w^ — —1, so that the specific enthalpy can be expressed in the following 
form 



Mm = yj-g ^^ • (21) 

One may also notice that 

Mm = W a $, a . (22) 

It is important to notice that the Clebsch potential $ has no direct physical meaning as it can be changed to another 
value $ — > $' = $ + $ such that the gauge function, $, is constant along the worldlines of the fluid: w a <b^ a = 0. 
In terms of the Clebsch potential the Lagrangian (19) of the ideal fluid is 



(23) 



Metrical tensor of energy-momentum of the ideal fluid is obtained by taking a variational derivative of the Lagrangian 
(23) with respect to the metric tensor, 

TSb - ■ (24) 



Calculation yields 

T aP = ( £ m + Pm) W a Wj3 + p m g a f3 , (25) 

where w a — dx a /dr is the four- velocity of the fluid, and r is the proper time of the fluid element taken along its 
worldline. This is a standard form of the tensor of energy- momentum of the ideal fluid [79]. Because the Lagrangian 
(23) is expressed in terms of the dynamical variable $, the Noether approach based on taking the variational derivative 
of the Lagrangian with respect to the field variable, can be applied to derive the canonical tensor of the energy- 
momentum of the ideal fluid. This calculation has been done in [58, pp. 334-335 ] and it leads to the expression 
(25). It could be expected because we assumed that the the ideal fluid consists of bosons. The metrical and canonical 
tensors of energy-momentum for the liquid differ, if and only if, the liquid's particles are fermions (see [58, pp. 331-332] 
for more detail). We do not consider the fermionic liquids in the present paper. 
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C. The Lagrangian of the Scalar Field 



The Lagrangian of the scalar field W is given by 

(26) 



where W = W is a potential of the scalar field. We assume that there is no direct coupling between the scalar field 
and the matter of the ideal fluid. They can interact only through the gravitational field. Many different potentials of 
the scalar field are used in cosmology [2]. At this step, we do not chose a specific form of the potential which will be 
selected later. 

Metrical tensor of energy-momentum of the scalar field is obtained by taking a variational derivative 



2 &C q 

that yields 



= (27) 



= d a Wp* - g a 



(28) 



The canonical tensor of energy-momentum of the scalar field is obtained by applying the Neother theorem and leads 
to the same expression (28). 

One can formally reduce the tensor (28) to the form similar to that of the ideal fluid by making use of the following 
procedure. First, we define the analogue of the specific enthalpy of the scalar field "fluid" 

Mq - V-g^d^d^ , (29) 
and the effective four- velocity, v a , of the "fluid" 

Mq«a = ■ (30) 

The four- velocity v a is normalized to v a v a = —1. Therefore, the scalar field enthalpy /i q can be expressed in terms 
of the partial derivative from the scalar field 

Mq = v a d a V . (31) 
Then, we introduce the analogue of the rest mass density p q of the scalar field "fluid" by defining, 



Pq EE ^ q = V a d a V = yZ-g^drVdv* . (32) 

As a consequence of the above definitions, the energy density, e q and pressure p q of the scalar field "fluid" can be 
introduced as follows 

e q = —g^d^d^ + Wm = i PqMq + W(^) , (33) 

Pq ee —g^d^d^ - Wm = i PqMq - Wm . (34) 



One notices that a relationship 



Mq = , (35) 

Pq 



between the specific enthalpy /z q , the density p q , the pressure p q and the energy density e q , of the scalar field "fluid" 
formally holds on the same form (12) as in the case of the barotropic ideal fluid. 

After applying the above-given definitions in equation (28), it is formally reduced to the tensor of energy- momentum 
of an ideal fluid 

T a/3 = ( £ q + Pq) U " W /3 + Pq9a/3 ■ (36) 

It is worth emphasizing that the analogy between the tensor of energy- momentum (36) of the scalar field "fluid" with 
that of the barotropic ideal fluid (25) is rather formal since the scalar field, in the most general case, does not satisfy 
all required thermodynamic equations because of the presence of the potential W = in the energy density e q , 

and pressure p q of the scalar field. 



12 



D. The Lagrangian of the Localized Astronomical System 



The Lagrangian £ p of matter of the localized astronomical system which perturbs the geometry of the background 
manifold of the FLRW universe, can be chosen arbitrary. We shall call the perturbation of the background manifold 
that is induced by £ p , the bare perturbation. We assume that the matter of the bare perturbation is described by a set 
of scalar potentials 9 which are analogues of the Clebsch potential of the matter supporting the background geometry. 
The Lagrangian density of the bare perturbation is given by L p = ^J—glP {6,g a @)- Tensor of energy-momentum of 
the matter of the bare perturbation, T Q( g, is obtained by taking a variational derivative 

2 



Tensor T Q/ 3 is a source of the bare gravitational perturbation of the background manifold that will be determined by 
solving Einstein's field equations derived in next sections. 



V. BACKGROUND MANIFOLD 



A. The Hubble Flow 



We shall consider the background universe as described by the Friedmann-Lemitre-Robertson- Walker (FLRW) 
metric. The functional form of the metric depends on the coordinates introduced on the manifold. Because the 
FLRW metric describes homogeneous and isotropic spacetime there is a preferred class of coordinates which clearly 
reveal these properties of the background manifold. These coordinates materialize a special set of freely falling 
observers, called comoving observers. These observers are following with the flow of the expanding universe and have 
constant values of spatial coordinates. The proper distance between the comoving observers increases in proportion 
to the scale factor R(T). In the preferred cosmological coordinates, the time coordinate of the FLRW metric is just 
the proper time as measured by the comoving observers. A particle moving relative to the local comoving observers 
has a peculiar velocity with respect to the Hubble flow. An observer with a non-zero peculiar velocity does not see 
the universe as isotropic. 

For example, the peculiar velocity of the solar system implies the dipole anisotropy of cosmic microwave background 
(CMBR) radiation corresponding to \vq\ — 369.0 ±0.9 km-s -1 , towards a point with the galactic coordinates (l,b) = 
(264°, 48°) [47, 53]. Such a solar system's velocity implies a velocity |t? LG | = 627±22 km-s" 1 toward (I, b) = (276°, 30°) 
for our Galaxy and the Local Group of galaxies relative to the CMBR [35, 55]. The existence of the preferred frame in 
cosmology should not be understood as a violation of the Einstein principle of relativity. Indeed, any coordinate chart 
can be used in order to describe the FLRW universe. A preferred frame exists merely because the FLRW metric admits 
only six-parametric group (3 spatial translations and 3 spatial rotations) as contrasted with the ten-parametric group 
of Minkowski (or De Sitter) spacetime which includes the time translation and three Lorentz boosts as well. The metric 
of FLRW does not remain invariant with respect to the time translation and the Lorentz transformations because 
its expansion makes different spacelike hypersurfaces non-equivalent. It may lead to some interesting observational 
predictions of cosmological effects within the solar system [60] . 



B. The Friedmann-Lemitre-Robertson- Walker metric 



In what follows, we shall consider the problem of calculation of the post-Newtonian perturbations in the expanding 
universe described by the FLRW class of models. The FLRW metric is an exact solution of Einstein's field equations 
of general relativity that describes a homogeneous, isotropically expanding or contracting universe. The general form 
of the metric follows from the geometric properties of homogeneity and isotropy of the manifold [99, 100]. Einstein's 
equations are only needed to derive the scale factor of the universe as a function of time. 

The most general form of the FLRW metric is given by 



ds 2 



-dT 2 + R 2 



dp 2 



l-kp 2 



p 2 (d 2 i? + sin 2 M 2 v) 



(38) 



where T is the coordinate time, {p, i9, v} are spherical coordinates, R = R(T) is the scale factor depending on time 
and characterizing the size of the universe compared to the present value of R = 1 . The time T has a physical meaning 
of the proper time of a comoving observer that is being at rest with respect to the cosmological frame of reference. 
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The present epoch corresponds to the value of the time T = Tq. The constant k can take on three different values 
k = { — 1,0, +1}, where k = — 1 corresponds to the spatial hyperbolic geometry, k = does the spatially flat FLRW 
model, and k = +1 does the spatially closed world [79]. 

The Hubble parameter H characterizes the rate of the temporal evolution of the universe. It is defined by 

IT.* = ±**. (39) 
R RdT y 1 

For mathematical reasons, it is convenient to introduce a conformal time, rj. via differential equation 

*-m- m 

If the time dependence of the scale factor is known, the equation (40) can be solved, thus, yielding T = T{rj). It 
allows us to re-express the scale factor R{T) in terms of the conformal time, R (T(rj)) = a(r)). The conformal Hubble 
parameter is, then, defined as 

:k^ = I* (4i) 

a adrj 

The two expressions for the Hubble parameters are related by means of equation 

H=^, (42) 

that allows us to link their time derivatives 

a 2 H = JC' - Vi 2 , (43) 
a 3 H = "K" - WK' + 2J{ 3 , (44) 

and so on. 

It is also convenient to introduce the isotropic Cartesian coordinates X 1 = {X,Y,Z} 7 by transforming the radial 
coordinate 

(45) 



k n 

l + -r 2 

4 



and defining r 2 = X 2 + Y 2 + Z 2 = 8ijX % X^ . In the isotropic coordinates the interval (38) takes on the following form 

ds 2 = G aP dX a dX 13 , (46) 
where the coordinates X a = {X° 7 X 1 , X 2 , X 3 } = {rj, X, Y, Z}, and the metric has a conformal form 

G aSi = a 2 (r])g a p (47) 



goo = -1 , goi = , gij 



I + -r 



(48) 



Determinant of the metric G Q( g is G = a 8 g, where g = — (l + fcr 2 /4) . The spacetime interval in the isotropic 
Cartesian coordinates reads 



ds 2 = a 2 (ri) 



, 2 8udX l dX^ 
-dr\ + ; 

k n 

1 + -r 2 

4 



(49) 



The distinctive property of the isotropic coordinates in the FLRW universe is that the radial coordinate r is defined 
in such a way that the three-dimensional space looks exactly Euclidean and null cones appear in it as round spheres 
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irrespectively of the value of the space curvature k. The isotropic coordinates do not represent proper distances on 
the sphere, nor does the radial coordinate r represents a proper radial distance measured with the help of radar 
astronomy technique. The proper spatial distance in the isotropic coordinates is (1 + kr 2 / '4)~ 1 ar [99]. 

The FLRW metric presented in the conformal form by equation (49) singles out a preferred cosmological reference 
frame defined by the congruence of worldlines of the fiducial test particles being at rest with respect to the spatial 
coordinates X 1 . Four-velocity of a fiducial particle is denoted as U a ~ dX a /dr, where dr = —ds is the proper time 
on the worldline of the particle. In the isotropic conformal coordinates, U a — (1/a, 0, 0, 0). The four-velocity is a unit 
vector, U a U a — G a ptJ a U 13 = —1. It implies that the covariant components of the four-velocity are U a — (—a, 0, 0, 0). 
In the preferred frame the universe looks homogeneous and isotropic. The choice of the isotropic Cartesian coordinates 
reflects these fundamental properties explicitly in the symmetric form of the metric (47). However, the set of the 
fiducial particles is a mathematical idealization. In reality, any isolated astronomical systems (galaxy, binary star, the 
solar system, etc.) have a peculiar velocity with respect to the preferred cosmological frame formed by the Hubble 
flow. We have to introduce a locally-inertial coordinate chart which is associated with the isolated system and moves 
along with it. Transformation from the preferred cosmological frame to the local chart must include the Lorentz boost 
and a geometric part due to the expansion and curvature of cosmological spacetime. It can take on multiple forms 
which originate from certain geometric and/or experimental requirements [16, 18, 48, 54]. 

We do not impose specific limitations on the choice of coordinates on the background manifold and keep the overall 
formalism of the post-Newtonian approximations, covariant. The arbitrary coordinates are denoted as x a — (x°,x 4 ) 
and they are related to the preferred isotropic coordinates X a = (rj^X 1 ) by the coordinate transformation x a = 
x a (X 13 ). This transformation has inverse X a = X a (x 13 ), at least in some local domain of the background manifold. 
In this domain, the matrices of the coordinate transformations 

A« 9 = JL, , m-„ . y , (50) 

and they satisfy to the apparent equalities A" 7 M 7 ^ = 5^ and M Q 7 A 7 ( g = 6p. 

Four-velocity of the Hubble observers written in the arbitrary coordinates has the following form 

u a = A Q (3 [7' 3 = a _1 A a , u a = MPjJp = -aM° a . (51) 

The background FLRW metric written down in the arbitrary coordinates, x a , takes on the following form 

g a p(x a ) = a 2 ± a/3 (x a ) . (52) 

Here the scalar function a(x a ) = a [r/(x a )], and the conformal metric 

f afi (x a ) = M^M^g^pr*) . (53) 

Any metric admits 3+1 decomposition with respect to a congruence of a timelike vector field [79]. FLRW universe 
admits a privileged congruence formed by the four- velocity u a of the Hubble observers which is a physically privileged 
vector field. The 3+1 decomposition of the FLRW metric is applied in arbitrary coordinates and has the following 
form 

9a/3 = -U a U/3 + P a f3 , (54) 

where the tensor 

P af} = a 2 M\W [iElJ , (55) 

describes the metric on the spacelike hypersurface being everywhere orthogonal to the four-velocity u a of the Hubble 
flow. Tensor P a p is the operator of projection on this hypersuface. It can be also interpreted as a metric on the 
hypersurace of orthogonality to the Hubble vector flow. Equation (54) can be used in order to prove that P a f) satisfies 
the following relationship 

pPvpp" = ; ( 56 ) 

which can be confirmed by inspection. The trace P a a — g a ^ ' P a p = P a ^P a p = 3. 

Now, we consider how to express the partial derivatives of any scalar function F = F(rj), which depends only on 
the conformal time r\ = rj(x a ), in terms of the four-velocity u a of the Hubble flow. Taking into account that rj = x° 
and applying equation (51), we obtain 

In particular, the partial derivative from the scale factor, a jCt = —au a = — !Ku Q , and the partial derivative from the 
Hubble parameter CK a = — 3~Cu a . 
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C. The Christoffel Symbols and Covariant Derivatives 

In the following sections of the paper we will need to calculate the covariant derivatives from various geometric 
objects on the background cosmological manifold covered by an arbitrary coordinate chart x a = (x°,x l ). The 
calculation engages the affine connection f ^ 7 of the background manifold which is decomposed into an algebraic 
sum of two connections (the Christoffel symbols) because of the conformal structure of the FLRW metric [98]. By 
definition, 

f a p~, = -j9 aiJ {g v p tl + g v ~t,i3 - 9mm) , (58) 

where 

= -2Hg a pu 7 + a 2 f Q/ 3 i7 . (59) 

Separating terms in the right side of (58) yields 

f a M = A a Pl + B a M , (60) 

where 

A a fil = -H{8$u 1 + 5«u fi - u a g Pl ) , (61) 

and 

= M/3>7 + ^7,/3 ~ * 07.J*) • ( 62 ) 

The non-vanishing components of the connections are given in the isotropic Cartesian coordinates X a by 

-o ~i k SpX q + S q X p - 5 pq X l 
A a ofi = JtSp , A ij = Wgij , B l pq = —- £ , (63) 

1 + -r 2 

4 

where X q = 5 q jX 3 , and all the other components of the connections vanish. 

A covariant derivative of a geometric object (scalar, vector, etc.) on the background manifold is denoted in this 
paper with a vertical bar. For example, the covariant derivative of a vector field F a is 

F a ]/s = F a ^ + T a 01 F\ (64) 

where a comma in front of sub- index (5 denotes a partial derivative with respect to coordinate . Equation (64) can 
be brought to yet another form if we denote the covariant derivative of the affine connection B a l g 7 with a semicolon. 
Making use of (60) in equation (64)transforms it to the following form 

F a ]/3 = F a , + A a ^ . (65) 

The covariant derivative of a covector F a is defined in a similar way, 

F a \p = F a ^ - t^pF^ (66) 

which is equivalent to 

F a \p = F a , - A^apFj , (67a) 

F a:JS = F a , p -ffT afi F y (67b) 

Equations for tensors of higher rank can be presented in a similar way. Of course, the covariant derivative of a scalar 
field F always coincides with its covariant derivative by definition, 

F\ a = F. a = F a . (68) 

We also provide an equation for the covariant derivative of the four- velocity of the Hubble flow. Doing calculations 
in the isotropic coordinates X a for the four- velocity U a , and applying the tensor law of transformation to arbitrary 
coordinates x a , results in 

u alfj = HP aP , u a ^ = H(6 a p + u a up) , u a ^ = HP"** , (69) 
where the tensor indices are raised and lowered with the metric g a p. 
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D. The Riemann Tensor 



The Riemann tensor is denned by 

R a i3fiv = ^ a pv,n — r° ji^. v + r^^r 7 ^ — r Q 1/7 r 7 ( 3 AJ . (70) 

and can be calculated directly from this equation. We prefer a slightly different way by making use of the algebraic 
decomposition of the Riemann tensor into the irreducible parts 

1 R 

Rafinv — Caji^y + — [S aix gp v + S[j„g a ^ — S^gp^ — S 0ll g av ) + — {g a ^(j[3v — (javgpn) , (71) 



where C a j3^ v is the Weyl tensor, 



S^v — R^v — -^Rg^v , (72) 



Rpv = g a ^ Raft/iv is the Ricci tensor, and R — g a ^R a /3 is the Ricci scalar. FLRW cosmological metric (38) has a 
remarkable property - it can be always brought up to the conformally-flat form by applying an appropriate coordinate 
transformation [49]. However, the Weyl tensor of any conformally-flat spacetime vanishes identically, 

Cafitiv = . (73) 

Direct evaluation of other tensors entering (71) by making use of the FLRW metric (47), (48) yields 

R„v - -^[Oi' (g^-2uf,u u ) + 2(^ 2 + k)(g^ + Uf,u u )] , (74) 
R = [JC + IK 2 + k] . (75) 
Making use of equations (73) - (75) in the decomposition (71) of the Riemann tensor, yields the following result 

R-afS^u = ^ {9an9Pv ~ 9a V 90n) - (W - 'K 2 - k) [P afl Pp„ - P av Pfi^)] , (76) 

where P a p = g a fj + u a up is the operator of projection that was introduced in (55). 



E. The Friedmann Equations 

The Einstein tensor £ a ^ = R a p — g a pR/2 of the FLRW cosmological model is derived from equations (74) and 
(75). It reads 

l a[i = - ^ [2 - K 2 - k) P a p + 3 (IK 2 + k) g a0 ] . (77) 
Einstein's field equations on the background spacetime takes on the following form 

l a p = 8,nT a p , (78) 

where the tensor of energy-momentum of the background spacetime manifold includes the background matter and 
the scalar field 

f a p=f™ +f* . (79) 

Here, tensors of energy-momentum in the right side of Einstein's equations are derived from the Lagrangians (23) 
and (26), and represent an algebraic sum of tensors (25) and (29). Each tensor of energy- momentum, T™ and T^p, 
is Lie-invariant with respect to the group of symmetry of the background FLRW metric independently, and each of 
them have the form of the tensor of energy-momentum of the perfect fluid. Hence, the tensor of energy-momentum 
T a fj in the right side of (78) has the form of a perfect fluid as well, 



T a fs = {e+p) u a u [3 + p g a/3 



(80) 
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It imposes a certain restriction on the effective energy density e and pressure p which must obey Dalton's law for 
a partial energy density and pressure of the background matter and the scalar field components [9] 

e = e m + e q , (81) 

P = Pm+Pq- (82) 

Here, e m and p m are the energy density and pressure of the ideal fluid, and e q and p q are the energy density and 

pressure of the scalar field which are related to the time derivative "J of the scalar field and its potential W = W(^>) 
by equations (33), (34). On the background spacetime these equations takes on the following form 

e q = \p q P q + W, (83) 

= ^PqMq- W, (84) 

where p q is the background specific enthalpy of the scalar field defined by (29), and p q = p q is the background density 
of the scalar field "fluid". It is worthwhile to remind to the reader that due to the homogeneity and isotropy of the 
FLRW universe, all matter variables on the background manifold are functions of the conformal time 77 only when 
being expressed in the isotropic Cartesian coordinates. 

Einstein's equations (78) can be projected on the direction of the background four-velocity of matter and on the 
spatial hypersurface being orthogonal to it. It yields two Friedmann equations for the evolution of the scale factor a, 



H>= ^e-*, (85) 
3 a z 



2H + 3H 2 = -87rp-4 (86) 
a z 

where e and p are the effective energy density and pressure of the mixture of matter and scalar field as defined above. 
A consequence of the Friedmann equations (85), (86) is an equation 

H- A = - 47 r(e+p) , (87) 

relating the time derivative of the Hubble parameter with the sum of the overall energy density and pressure, which 
can be expressed in terms of the density and specific enthalpy of the background components of matter, 

£ + P = PmMm + PqPq ■ (88) 

In order to solve the Friedmann equations (85), (86) we have to employ the equation of state of matter. Customarily, 
it is assumed that each matter component obeys its own cosmological equation of state, 

Pm = w m e m , p q = w q e q , (89) 

where w m and w q are parameters lying in the range from —1 to +1. In the most simple cosmological models, 
parameters w m and w q are fixed. More realistic models admit that the parameters of the equation of state may 
change in the course of the cosmological expansion, that is they may depend on time. The equation of state does not 
close the system of the Friedmann equations, which have to be complemented with the equations of motion of the 
scalar field and of the ideal fluid in order to make the system of differential equations for the gravitational and matter 
field variables complete. 



F. The Hydrodynamic Equations of the Ideal Fluid 

The background value of the Clebsch potential of the ideal fluid, $, depends only on the conformal time ij of the 
FLRW metric. The partial derivative of the potential, taken in arbitrary coordinate chart on the background manifold, 
can be expressed in accordance with equation (57) in terms of the background four-velocity u a as follows 

#| Q = -$u a . (90) 

It allows us to write down the specific enthalpy of the ideal fluid in terms of the Clebsch potential. Taking background 
value of equation (22), we obtain 



(91) 
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The background equation of continuity for the rest mass density p m of the ideal fluid is 

(p m u a ) [a = , (92) 

that is equivalent to 

p m \ a - 3Hp m u a = , (93) 
where we have used (69). The background equation of conservation of energy is 

e m |a ~ 3# (e m + p m ) u a = , (94) 
where we have employed definition of the energy (11), and equation (93) along with (13). 

G. The Scalar Field Equations 

Background equation for the scalar field # is derived from the action (8) by taking variational derivatives with 
respect to ty. It yields 

a- dW 

?**\<# -qW = ° • 

In terms of the time derivatives with respect to the conformal time 77, equation (95) reads 

dW 

* + 3ff * + — — = . (96) 

Here, we have taken into account that the background value of the scalar field, depends only on time 77, and its 
derivative (with respect to 77) is proportional to the background four-velocity 

*|a = -*«a • (97) 

If we use definition of the background enthalpy of the scalar field 

,I q = u a $ la = 4> , (98) 

and account for definition (33) of the specific energy e q of the scalar field, the equation (96) will become 

e qla ~3H{e q +p q )u a = . (99) 

that looks similar to the hydrodynamic equation (93) of conservation of energy of the ideal fluid. Because of this 
similarity, the second Friedmann equation (86) can be derived from the first Friedmann equation (85) by taking a 
time derivative and applying the energy conservation equations (94) and (99). 

The background density p q of the scalar filed "fluid" is p q = pZ q in accordance with (32). The equation of continuity 
for the density p q of the ideal fluid is obtained by differentiating definition of p q , and making use of (96). It yields 

dW 

fo fia )i« = -w (100) 



or, equivalently, 



dW _ 

p q \ a - 3Hp q u a = -g^-Ua , (101) 



which shows that the "density" of the scalar field "fluid" is not conserved. We emphasize that there is no any 
violation of physical laws, since (101) is simply another way of writing equation (95), and the scalar field is not 
thermodynamically equivalent to the ideal fluid. Equation (101) is convenient in the calculations that follows in next 
sections. 
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H. Equations of Motion of Matter of the Localized Astronomical System 

Matter of the localized astronomical system is described by the tensor of energy-momentum 1 a P defined in (37) in 
terms of the Lagrangian derivative. It can be given explicitly as a function of field variables after we chose a specific 
form of matter, for example, gas, liquid, solid, or something else. We do not restrict ourselves with a particular form 
of this tensor, and shall develop a more generic approach that is applicable to any kind of matter comprising the 
localized astronomical system.. 

Background equation of motion of matter of the astronomical system is given by the conservation law 

i q V = 

It can be also written down in terms of a covariant derivative of the conformal metric 

where the connection A a J g 7 is defined in (61). 

It is natural to write down this equation in 3+1 form by projecting it on the direction of 4-velocity of the Hubble 
flow, u a , and on the hypersurface being orthogonal to it. This is achieved by introducing the following projections 



a = u^u u 1^ , (104a) 

r = P^V , (104b) 

r a = -P^u'l^ , (104c) 

Ta(i = P^Pf^l^, (104d) 



which corresponds to the kinemetric decomposition of 1^ introduced by Zelmanov [108]. Quantity a is the energy 
density of matter of the localized system, t a is a density of linear momentum of the matter, and t a p is the stress 
tensor of the matter. 

Equations of motion (102) of the localized matter can be rewritten in terms of the kinemetric quantities as follows, 

Ht , (105a) 
H (t q - u a T) . (105b) 

Equation (105a) is equivalent to the law of conservation of energy of matter of the localized system. Equation (105b) 
is analogues to the Euler equation of motion of fluid or the equation of the force balance in case of solids. 

VI. LAGRANGIAN PERTURBATIONS OF FLRW MANIFOLD 

A. The Concept of Perturbations 

In the present paper, FLRW background manifold is defined by the metric g a p which dynamics is governed by 
background matter fields - the Clebsch potential <l of the ideal fluid and the scalar field We assume that the 
background metric and the background values of the fields are perturbed by a localized astronomical system which 
is considered as a bare perturbation associated with a field variable 0. Perturbations of the metric and the matter 
fields caused by the bare perturbation are considered to be small so that the perturbed metric and the matter fields 
can be split in their backgrounds values and the corresponding perturbations, 

g a p = g a p + x a p , $ = $ + 4> , # = * + ip ■ (106) 

These equations are exact. We emphasize that all functions entering equation (106) are taken at one and the same 
point of the background manifold. The bare perturbation does not remain the same in the presence of the perturbations 
of the metric and the matter fields. Therefore, the field variable corresponding to the bare perturbation, is also 
perturbed 

6 = 6 + 6 . (107) 

We consider the perturbations of the metric - x a p, the Clebsch potential - (f>, and the scalar field - ip as being weak 
with respect to their corresponding background values g a p, and 'J, which dynamics is governed by the background 



(102) 



(103) 



(<™ a +r«) |Q = - 
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equations that have been explained in section V. Because the field variable 9 is the source of the bare perturbation, 
we postulate that its background value is equal to zero: = 0. The perturbations >t a fs, <j>, and ip have the same order 
of magnitude as 9. 

Perturbation of the contravariant component of the metric is determined from the condition g ai g 713 = g ai g lfi = S^, 
and is given by 

g afj = g afi - x al) + ^ Q 7 x^ + . . . , (108) 

where the ellipses denote terms of the higher order. 

It turns out [44, 85] that a more convenient field variable of the gravitational field in the the theory of Lagrangian 
perturbations of curved manifolds, is a contravariant (Gothic) metric 

Q ap = V~99 aP ■ (109) 

The convenience of the Gothic metric stems from the fact that it enters the de Donder (harmonic) gauge conditions 
which significantly simplifies the Einstein equations [65, 98]. The Gothic metric variable is also indispensable for 
concise and elegant formulation of dynamic field theories on curved manifolds [26]. Making use of the Gothic metric 
allows us to significantly reduce the amount of algebra in taking the first and second variational derivatives from the 
Hilbert Lagrangian and the Lagrangian of the background matter in FLRW universe as explains in the rest of this 
section. 

The covariant Gothic metric is defined by means of equation 

a V 7 = 6° , (110) 
that yields g Q( g = g a fj/\/—g- We accept that Q a ^ is expanded around its background value, Q a @ — ^/—gg al3 , as follows 

which is an exact equation. 

Further calculations prompt that it is more suitable to single out ^J~-g from t)"' 3 , and operate with a variable 

i°p = l_ . (n2) 

This variable splits the dynamic degrees of freedom of the gravitational perturbations from the background manifold 
which evolves in according with the unperturbed Friedmann equations. Tensor indices of are raised and lowered 
with the help of the background metric, for example, l a p = g a ^g0v^ v ■ The field variable l a P relates to the perturbation 
K a p of the metric tensor. To establish this relationship, we start from (109), substitute equation (111) to its left side, 
and expand its right side in the Taylor series with respect to x a p. It results in 

"" f = ^^ + 2W^'' + ---- (U3) 
where the partial derivatives are calculated by successive application of the following rules 

= - V f^{g°" i Sr Bv + 9°"^ -9 aP aT) , (H4a) 
f)a a P 1 

J— = -^(r'l 1 "" , (iwb) 

W = • (n4c) 

which can be easily confirmed by inspection. Replacing the partial derivatives in (113) and making use of the definition 
(112), yields the relationship between l a/B and x a/3 as follows 

l"P = ~x af3 + \g a(i K + yt^ a ^\ - ! x - ^g af3 (x^x^ - \^ + ... , (115) 

where h = „ = g pa *c p<7 , and ellipses denote terms of the cubic and higher order in x a p. 
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Perturbations of four- velocities, w a and v a , entering definitions of the energy- momentum tensors (25), (36), are 
fully determined by the perturbations of the metric and the potentials of the matter fields. Indeed, according to 
definitions (20) and (32) the four-velocities are defined by the following equations 

w a = , v a = . (116) 

Mm Mq 



where fi m — y / ~g a i 3 ^^jj and ^ q = — g ali ^ ^ ,p in accordance with (21) and (29) respectively. We define 
perturbation of the covariant components of the four- velocities as follows 

W a — U a + Sw a , V a =Ua+Sv a , (117) 

where the unperturbed values of the four- velocities coincide and are equal to the four- velocity of the Hubble flow due 
to the requirement of the homogeneity and isotropy of the background FLRW metric. Substituting these expansions 
to the left side of definitions (116), and expanding its right side by making use of the expansions (106) and (108) of 
the scalar fields and the metric, yields 

Sw a = -^-P ? a ^\/3 - \(\u a , 5v a = -^P P a ip\p - \<\u a , (118) 

Mm ^ Mq ^ 



where we have introduced a new notation 



q = -u a u p x a p , (119) 



for the gravitational perturbation of the metric tensor projected on the background four- velocity of the Hubble flow. 
Making use of l a p , the previous equation can be recast to 

q = u a uH afi + l - , (120) 

where I = l a a — g af) l a fj. Remembering that g a/3 = P Q/3 — u Q u /3 , we can put equation (120) yet to another form 

q = \ (u a u P + P aP ) l a p , (121) 
which is useful in the calculations that follows. 



B. The Perturbative Expansion of the Lagrangian 



We have introduced the Lagrangian of the theory in section IV. The Hilbert Lagrangian of the gravitational field 
is £ s = — ^/^gR/16ir, where R is the Ricci scalar. The Lagrangian density of matter is L m = y/^g~L m ($>,g a p), and 
the Lagrangian density of the scalar field L q = ^/^gL q (^, g a p). The matter, the scalar field as well as the spacetime 
manifold are perturbed by a matter of N-body system described by a set of field variables O with the Lagrangian 
density L p = ^gL v (Q, g a p). 

The action of the unperturbed FLRW universe is a functional 

S = / d A xl , (122) 
Jm 

depending on the unperturbed Lagrangian 

£=Z s +£ m +£ q , (123) 

taken on the background values of the field variables g a p, #, and 

The presence of a localized astronomical system perturbs the spacetime manifold and the background values of the 
field variables. The perturbed Lagrangian becomes an algebraic sum of four terms 

L =L S +L m + L q + L p , (124) 

where the Lagrangian £ p describes the bare perturbation, and L s , £ m , L q are perturbed values of the Lagrangian of 
the FLRW universe. 
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The perturbed Lagrangian can be decomposed in a Taylor series with respect to the perturbed values of the field 
variables. It is achieved by substituting expansions (106) to the Lagrangian (124) and expanding it around the 
background values of the variables. It yields 

L = L+L 1 +L 2 +L 3 + ... , (125) 

where L\, L 2 , L 3 , ■ . ■ , are the Lagrangian perturbations which are linear, quadratic, cubic, and so on, with respect 
to the perturbations of the field variables, h a p, <p, ip, and 9. More specifically [85], 

Ll = ^^_ + ^^+^" +(C p ) (126a) 

- ^HIK4(4K^(*f) ,126b) 



and so on. Here, the variational derivatives from the Lagrangian density, L, depending on the field variables and their 
derivatives, are defined as follows 

d 2 dL 

(127a) 



SL 


dL 


d dL 




d^ v 


dx a dg>% 


SL 


dL 


d dL 


<5<l 


9$ 


dx a <9<& Q ' « 


SL 


dL 


d dL 


~sH 




dx a d^. a 



dx a dx$ dg^ iCe 
d 2 dL 



dx a dxP d$, a0 ' 

d 2 dL 
dx a dxf } d^ ~r ' 



(127b) 
(127c) 



The variational derivative with respect to the metric density g^ v relates to the derivative with respect to the metric 
by an algebraic operator 

One has to notice that the expansion (125) is defined up to the terms which are represented as a total covariant 
derivative from a vector density (the, so-called, divergent terms). For example, the direct Taylor series expansion 
shows that the Lagrangian £ 2 has a term with cross-coupling of (f> and i[>. This term was eliminated from £ 2 because 
it can be represented as a covariant divergence from a vector that vanishes identically after taking the Lagrangian 
derivative [85]. The divergency terms can be important in the discussion of the boundary conditions but they do not 
enter the equations of motion of fields which represent a system of differential equations in partial derivatives for the 
perturbations of the dynamic (field) variables. Furthermore, it is straightforward to prove that any of the Lagrangian 
derivatives (127a)-(127c), applied to a partial derivative of a geometric object F = F(Q a P; $; g al3 17 ; $ )7 ; ^ ;7 ; . . .), 
vanishes [80] 

dF ' 0, 'm-o, x(s=o. 



5Q a P \dx a J ' <S$ \dx a J ' 5^ \dx a 

Equations (129) does not hold for a covariant derivative [80]. We shall use equation (129) for bringing the Lagrangian 
derivatives to a simpler form. 

The field equations are obtained by taking the variational derivatives from the perturbed action with respect to 
various variables subject to the least action principle. In accordance with this principle, the variational derivatives 
from the perturbed Lagrangian must vanish, 

M ~ SJC n SL n 

= 0, — = , — =0 . (130) 



(5$ ' (5* 

We substitute the Taylor decomposition (125) of the Lagrangian to equations (130) and separate the background value 
of the derivatives from their perturbed values. We assume that gravitational dynamics the unperturbed universe obeys 
the background field equations. Then, the perturbed part of the equations represent a series of equations of the first, 
second, third, etc. order, which can be solved by successive iterations. In this paper we restrict ourselves with the 
linearized approximation of the first order with respect to the perturbations. It generalizes the first post-Minkowskian 
approximation to the case of the expanding universe. 
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C. The Background Field Equations 

The dynamics of the background universe is governed by the variational equations 



SL g SL m SL q n , 
+ + = °> ( 131a ) 

O0Q/3 



= , (131b) 
= . (131c) 



After performing the derivatives, equation (131a) becomes the Einstein equation (78), equation (131b) is reduced to 
equation of continuity (92) after taking into account the thermodynamic relationship (15), and equation (131c) is 
equivalent to (95). These equations have been thoroughly discussed in section V. Solution of these equations depend 
on equation of state of background matter. We assume that the solution exists and that the time dependence of the 
FLRW metric g a p — g a /3(v)> the Clebsch potential $ = $(??), and the scalar field ^> = ^(r/) is explicitly known. 



D. The Gravitational Field Perturbations 

The gravitational field perturbations satisfy the following (exact) differential equation 

F„„ = 8^ ( V + V) . (132) 

which generalizes the Einstein field equations of the post-Minkowskian approximations in asymptotically flat spacetime 
to the case of the expanding universe. Tensor is an algebraic superposition 

F^ = F^ + F^+F^, (133) 

where the linear operators in the right side are defined through the Lagrangian derivatives as follows, 



g bg^ v \ 5g 

F m = -i^J-lV^+M^ (134b) 
FS, .-J£, ' (134c, 



-gSg^ V h af3 

The right side of equation (132) contains the tensor of energy-momentum of the bare gravitational perturbation 
which is generated by the matter of the localized astronomical system and can be calculated as the Lagrangian 
derivative (37) . The right side of (132) also contains the non-linear corrections that are given by 

t^ = ^( — + —!- + ...) . (135) 



-g \Sg^ 8g~v-v 

In what follows, we shall neglect the contribution of t M „ as it is of the higher order compared with other terms in 
(132). 

The differential operator, Fj* v , represents a linearized perturbation of the Ricci tensor, and after calculation of 
(134a), is given by 

= 2 (jv^ \a 9n"l ^ \ctf3 ~ la[i\v" ~ ^av\n ^ ? (136) 

where each vertical bar denotes a covariant derivative with respect to the background metric g^. 

Operators i* 1 ™ and Fj} v depend essentially on a particular choice of the Lagrangian of matter and scalar field, and 
take on different forms depending on the specific analytic dependence of L m and L q on the field variables. In the 
particular case of the ideal fluid, the term embraced in the round parentheses in the right side of equation (134b) is 

f)Q/3 ^ + lT = ^lT%--g a0 f™)+cpd a (p in V=3u a ) , (137) 
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where u a = —g al3 <&^/p m , and 7™g is given in (25). We emphasize that though the ideal fluid satisfies the equa- 
tion of continuity (92), it should not be immediately implemented in (137) because this expression is to be further 
differentiated with respect to the metric tensor according to (134b). 

For the scalar field, the term enclosed to the round parentheses in the right side of (134c) is 



= 0\V _ ._ 



-gu 



(138) 



where u a = —g a ^^p/p q , p q — p q , T^a is given in (36), and the equation of continuity for the scalar field (100) should 
not be implemented until differentiation with respect to the metric tensor (134c) is completed. 

Taking the variational derivatives with respect to g^ v from the expressions (137) and (138), and applying thermo- 
dynamic equations (18), allows us to write down the right sides of equations (134b), (134c) as follows, 



F m = -4tt 



(p 



+8np„ 



1 



(e m +Pm)qu tl u lJ 



(139) 



1 



-47T 



(Pq 



dW , 

e q )l fny -2g^-^-ij 



+ 87rp q {u^ tV + - g^w Q V,a) , 



(140) 



where p q = ^ /a in accordance with definition (32). The potential energy of the scalar field, W = W(^>) 7 remains 
arbitrary as yet. 

It is important to emphasize that in the most general case the ratio Cg /c 2 of the speed of sound in fluid to the 
fundamental speed c, is not equal to the parameter w m of the equation of state (89), that is w m ^ (c s /c) 2 . Indeed, 
the speed of sound is defined as a partial derivative of pressure p m with respect to the energy density e m taken under 
the condition of a constant entropy s m , 



dp n 



This equation is equivalent to the following relationship 



s m — const. 



(dp m /dp m ) s 



(de m /dp m ) s 



(141) 



(142) 



which is a consequence of thermodynamic relationships and definition of a partial derivative. The ratio of the partial 
derivatives in (142) is not reduced to u> m in case when w m depends on some other thermodynamic parameters implicitly 
depending on the specific enthalpy. For example, in case of an ideal gas the equation of state p m = w m e m , where 
w m = kT/mc 2 , k is the Boltzmann constant, m - mass of a particle of the ideal fluid, and T is the fluid temperature. 
The speed of sound c 2 = c 2 (dp m /de m ) s =const = Tw m > w m = p m /e m , where T > 1 is the ratio of the heat capacities 
of the gas taken for the constant pressure and the constant volume respectively. 

The scalar field with the potential function W{^) ^ does not bear all thermodynamic properties of an ideal fluid. 
Nevertheless, we can formally define the speed of "sound" c s propagating in the scalar field "fluid", by equation being 
similar to (142). More specifically, 



@Pq/fyq)*=c 



c 2 {de q /d^ =coQst ' 

Simple calculation reveals that the speed of "sound" of the scalar field is always equal to the fundamental speed 



(143) 



c s c , 



(144) 



irrespectively of the value of the potential function W(^). It explains why the terms being proportional to the factor 
1 — c 2 /c 2 , do not appear in the expression (140) as contrasted with (139). 
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E. The Ideal Fluid Perturbations 

The perturbed field equations for the ideal fluid are obtained by taking the variational derivatives with respect 
to the field <& from the Lagrangian (124) - it corresponds to the middle equation in (130). Taking into account the 
background equation (131b) yields the equation of sound waves propagating in the fluid as small perturbations, 

F m = 87rS m , (145) 

where the linear differential operator 

^ S + (146) 



-g8<& \ Sq^ 

and the source term 



1 / 5Lf SL^ 

S m = - — = —J- + —J- + ... . (147) 



In the case of a single-component ideal fluid, the Lagrangian (23) depends merely on the derivative of the Clebsch 
potential <& and on the metric tensor. Therefore, the explicit form of the linear operator F m is reduced to a covariant 
divergence 

F™ = Y a [a , (148) 

where a vector field 



Y c 



<93> 



1 -~2 lg ) [d^-2 9 ^ L ) +4> ' fi W„ 



(149) 



where the partial derivatives are taken from the Lagrangian L m , but not from its density £ m = y/—gL m . More 
specifically, calculations yield 

Y ^ 0|c _ p m r/3 % + ( i _ SL\ (^u% p - ^ m ^q) . (150) 

Similar equations were derived by Lukash [73] who used the variational method to analyze production of sound waves 
in the early universe. 

F. The Scalar Field Perturbations 

Equations for the scalar field perturbations are derived by taking the variational derivative from the Lagrangian 
(124) with respect to the field variable VP - see the last equation in (130). Subtracting the background equation (131c) 
leads to 



where the linear differential operator 



and the source term 



F q = 87rS q , (151) 



^ = T^ ('# + # + -) • (153) 



According to equation (26), the Lagrangian density of the scalar field £ q = v / — 5-k q depends on both the field and 
its first derivative, $f a . For this reason, the differential operator F q is not reduced to the covariant derivative from a 
vector field as the partial derivative of the Lagrangian with respect to VP does not vanish. We have 

I dW d 2 W 
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where I = g a/3 l a /3, and vector field 



Z a = 



V"' 



Performing the partial derivatives in equation (155), yields a rather simple expression 



Z a = ^ a - p q l ap u p 



a/3- 



(155) 



(156) 



where we have used equation = —u"tyipu a — —p q u a . The reader is invited to compare equation (156) with 
(150) to observe the differences between the Lagrangian perturbations of the ideal fluid and the scalar field. One may 
observe that (150) becomes identical with (156) in the limit c s — > c, and p m — > £L m . This corresponds to the case of 
an extremely rigid equation of state w m — 1 in equation (89). According to the discussion following equations (143), 
(144) the speed of 'sound' c s in the scalar field 'fluid' is always equal to c. However, it does not assume that the 



parameter w q of the equation of state of the scalar field, p q = 



in (89) is equal to unity. This is because the 



scalar field is not completely equivalent to the ideal fluid in the sense of thermodynamic [2]. . 



G. The Lagrangian Equations for Field Variables 

1. Equations for the metric tensor perturbations 



Linearized equations for gravitational field variables, l^, are obtained from (131a) after neglecting in its right side 
the non-linear source t M „, and rendering a series of transformations which re-arrange and sort out similar terms. First, 
let us make use of Einstein's equations (139) and (140) to find 



(157) 



+ 8?rp m 
+ 87rp q 



Un<j>, v + u„(j) ttl - g^u a (f>, a + 1 - 



'.a - -^Pmq I U^U V 



u M Vv + u u ip itl - g ttv u a tl) ta + g^ 



9* /iq 



Second step is to transform the linear differential operator Ff> v in (136) to a more convenient form that will allow 
us to single out the gauge-dependent terms denoted by 



(158) 



Changing the order of the covariant derivatives in (136) and taking into account that the commutator of the second 
covariant derivatives is proportional to the Riemann tensor, we can recast (136) to the following form, 



fib 1 



1 



(159) 



where the round brackets around indices denote symmetrization. The terms with the Ricci and Riemann tensors can 
be expressed in terms of the total background energy and pressure of the ideal fluid and scalar field by making use of 
equations (71), (73) and Einstein's equations (78). It yields 



R a ( t ih)u + R f iai3vl at3 = 47T 



p) V + 7i [P~ |) 9nv + (e+p) (2(/"(/ / ,/, /n - •2u n ii„l /ln (/,,(/„/ - ii iu a]) 



Finally, substituting equations (157), (159) and (160) to (132) results in 



/ \ct \ T, A a A A 



(160) 



(161) 



-167T 



£ l_ 
3"" + 4 



+167rp n 



<!,,„ f (< ■- ! a" n :: l:.,. 4 , ^"/."•■ / 

1 



u v qb^ - g^u a (j) ia + 1 - 



+167rp q 



+ u u ip^ - g^„u a ^ ia + g^ v 



p q 
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where the non-linear term was neglected. 

The first term in (161) is a covariant Laplace-Beltrami operator, lfj,J a \ a = g a ^l ft ^\ a p, that is a rather complicated 
geometric object. Its explicit expression can be developed by making use of the Christoffel symbols given in (60). 
Tedious but straightforward calculation yields 



l fj,u \a 



V ; % + 2Hu a l^ ;a - 2 (Hu a l a J lv - 2 

+ 2E(u^A v + u v An) + 2H 1 (l^ u — u a u^l 

+ 2H 2 (21 + 3u^u a l av + 3u„u a l afl - g^u a u^l a0 - u^uj) , 



(162) 



where the semicolon denotes a covariant derivative that is calculated with the Christoffel symbols B a ^ like in (67b), 
and the differential operator i M „ ;a ;Q = g al3 l^- a p- 

Further derivation of the differential equations for linearized metric tensor perturbations can be significantly sim- 
plified if we re-define the gauge function, A a = in the following form 



A a = -2Hl al3 up + 16tt (p m <P + pjj)) u a + B c 



(163) 



where B a is an arbitrary gauge function. This choice of the gauge function A a allows us to eliminate two terms in 
equation (162) which depend on the first covariant derivatives with respect to the background metric g a p. Moreover, 
it allows to eliminate a number of terms depending on the first derivatives of the fields <j> and i/j in equation (161). 
Since we keep the gauge function B a arbitrary, the equation (163) does not fix any gauge. The choice of the gauge is 
controlled by the gauge function B a . 

One substitutes the gauge function (163) to equations (162) and (161) and make use of the background Friedmann 
equations (85), (86) to replace the background values of the energy density, e, and pressure, p, with the Hubble 
parameter H and its time derivative H ' . It brings about equation (161) to the following form 



l^' a . a + 2Hu a l^. a + 2(ir + H 2 ) (J„„ + u^u a l au + u v u a l 
2k r 

2 l\iv + 2U[iU lew + 2liu1l l a ^ 



-rT67ra M u„ 



-~2^ 



dW 



+gtivB a \a — Bv\i> — Bv\(x + 2H(u ll B u 4- u y B j 



4> + Pqlp) 

j, - gf, u u a B a ) 



(164) 
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This equation is fully covariant and is valid in any gauge and/or coordinate chart. It clarifies the advantage in the 
choice of the gauge function (163). 

Indeed, if one works in the isotropic coordinates associated with the Hubble flow, where U a — (1/a, 0, 0, 0), it allows 
us to fully decouple the differential equations for Z o, ^oi and kj components of the metric tensor perturbations. Let 
us assume, for simplicity, B a = that is an analogue of the harmonic gauge in asymptotically-flat spacetime. Then, 
different tensor components of equation (164) become 



Dq + 23iq.. + 4kq - 4tt ( 1 - J p m p m q = 8tt (T o + Ifcfc) 



8na 



Pn 



c 2 



dW 



ni 



<ij> 



2HO<jj>;0 - 



f 23-Cloi-.o + 2kloi 
2 (%' - k) l Kij> 



Ul kk + 2m kk , + 2 (%' + 2k) 



k k 



167tToj , 
167rT<ij> , 
167rT fe fc ■ 



(165a) 



(165b) 
(165c) 
(165d) 



Here, we denoted Dl^ = f^l^-ap, q = (loo + hk) /2, hk = ^11 + ^22 + ^33, l<ij> = Uj-Q/tySijlkk, and the same index 
notations are applied to the tensor of energy-momentum of the localized astronomical system. These equations 
are clearly decoupled from one another, thus, demonstrating the advantage of the gauge condition (163) used along 
with B a = 0. 

Equations (165b)-(165d) can be solved independently if the initial and boundary conditions are known, and the 
tensor of energy-momentum of the localized astronomical system is well-defined. Equation (165a) for scalar q besides 
knowledge of 1 a 0, demands to know the scalar field perturbations, <f) and i[>, that contribute to the source of q standing 
in the right side of (165a). Equations for these perturbations are obtained in the following text. 
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2. Equations for the ideal fluid perturbations 

The ideal fluid perturbations, </>, evolve in accordance with the Lagrangian equation (145). In the linear approxi- 
mation we can neglect the source term S m in its right side. The covariant derivative in the definition of the linear 
operator F m given by (148), can be explicitly performed that yields equation for the Clebsch potential 

- 2fl m Hq + 167r/2 m (p m </> + p q ip) + ^1 - ^ ^u a u p <t>\ a p - ^ m « a q, a j = p, m u a B a , (166) 

where the gauge function (163 has been used. The gauge B a remains yet unspecified so that equation (166) is covariant 
and is valid in any coordinate chart. 



3. Equations for the scalar field perturbations 

Linearized equation for the scalar field perturbations, ip, is obtained from the Lagrangian equation (151) after 
neglecting the (non-linear) source term S q . After performing the covariant differentiation in equation (154), we get 
equation for the scalar field perturbation 

/ _ dW\ _ _ d 2 W 

where equation (163) has been used along with the equality p q = p, q . The gauge function B a is kept unspecified so 
that equation (166) is covariant and is valid in any coordinates. 



VII. GAUGE-INVARIANT FIELD EQUATIONS IN 3+1 FORMALISM 
A. Algebraic Decomposition of the Metric Perturbations. 

We have derived the system of coupled differential equations (164), (166), (167) for the field variables l a p, <f> and tp, 
describing perturbations of the gravitational field, the ideal fluid, and the scalar field respectively. These system of 
equations can be split into a set of gauge-invariant differential equations for the scalar, vector, and tensor parts which 
is convenient for theoretical study of the evolution of the perturbations in arbitrary coordinates. This 3+1 split is 
achieved by making use of the operator of projection P a p onto a hypersurface being orthogonal to the congruence of 
worldlines of the Hubble flow. 

The theory under development admits four, algebraically-independent scalar perturbations. Two of them are the 
Clebsch potential of the ideal fluid <j> and the scalar field ip. The two other scalars characterize the scalar perturbations 
of the gravitational field. They can be chosen, for example, as a projection of the metric tensor perturbation on the 
direction of the background four- velocity field, u a u^l a p, and the trace of the metric tensor perturbation, I = g a/3 l a p. 
However, it is more convenient to work with two other scalars, defined as their linear combinations, 

q = l - {u a uP + P a ' 3 ) l aP , (168a) 

p = (u a u?+g a0 )l a , (168b) 

Notice that the scalar q has been introduced earlier in (140). The scalar p is, in fact, projection of l a p onto the space- 
like hypersurface being orthogonal everywhere to the worldlines of fiducial observers moving with the background 
four-velocity u a of the Hubble flow. Indeed, after accounting for definition (54), equation (168b) can be written as 

P = P af3 l a p ■ (169) 

Vectorial gravitational perturbations are defined by a spacial-temporal projection 

p a = -Pj&lfr , (170) 

where minus sign was taken for the sake of mathematical convenience. Due to its definition, vector p a — g a ^p/3 
is orthogonal to the four- velocity u a , that is u a p a = 0. Hence, it describes a space-like vector-like gravitational 
perturbations with three algebraically-independent components. 
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Tensorial gravitational perturbations are associated with the projection 

Plli = Pa? - \P*P$ , (171) 

where 

p af} = P^Pfl^ . (172) 

Here, the tensor p Q( g is a double projection of l a p onto space-like hypersurface being orthogonal to the worldlines of 
fiducial observers moving with the four- velocity u a of the Hubble flow. The trace of this tensor coincides with the 
scalar p. Indeed, 

3 Q V/3 = g^P^P^i^ = p^p^i^u = P^ipv = p , (173) 

where the property of the projection tensor P^^Pp" — P^ v has been used. Equation (173) makes it clear that tensor 
p 7 a p is traceless, that is g a/3 p J a p = 0. Because of this property, and four orthogonality conditions, u a p J afj = 0, the 
symmetric tensor p J a o has only five, algebraically-independent components. 

Gravitational perturbation l a p can be decomposed into the algebraically-irreducible scalar, vector and tensor parts 
as follows 

la/3 = Pip + U a pp + Upp a + (u a Up + T^Papj P + 2u a Uj3 (q - p) . (174) 

One should not confuse the pure algebraic decomposition of the metric tensor perturbation with its decomposition 
in a functional (Hilbert) space. This decomposition was pioneered by Lifshitz and Khalatnikov [70] and later on, 
structured by Arnowitt, Deser and Misner Arnowitt et al. [5] (see also [79]). It is commonly used in the research on 
the relativistic theory of formation of the large-scale structure in the universe. The functional ADM decomposition 
of the metric tensor perturbations is done with respect to the direction of propagation of weak gravitational waves 
and singles out the longitudinal (L), transversal (T) and transverse-traceless (TT) parts of the perturbations. In 
other words, the functional decomposition make sure that the vector p a and the tensor parts of the gravitational 
perturbation, p T a p, are further decomposed in the functionally- irreducible components that are reduced to two more 
scalars, and two transverse vectors each of which has only two, functionally-independent components. The remaining 
part of the tensor perturbations, p T a a, is transverse-trackless and has only two functionally- independent components 
denoted as p^. The ADM decomposition of the metric tensor is a powerful technique in the theory of gauge-invariant 
cosmological perturbations [8] . However, it is not convenient in the development of the post- Newtonian dynamics of 
celestial bodies in cosmology, and shall not be used in the present paper. 

Our next step is derivation of the field equations for the algebraically-irreducible components of matter and gravita- 
tional field. Before doing this derivation, let us discuss the gauge transformations of the corresponding field variables. 



B. The Gauge Transformation of the Field Variables 

Gauge invariance is a cornerstone of the modern theoretical physics with a long and interesting history [52]. Gauge 
invariance should be distinguished from the coordinate invariance or the general covariance because, by definition, a 
gauge transformation changes only field variables of the theory under consideration but not coordinates. Discussion 
of gauge transformation and invariance requires introducing a gauge field and a new geometric object - an affinc 
connection - on a fibre bundle manifold describing the intrinsic degrees of freedom of corresponding field variables of 
the gauge field theory. 

The present paper discusses physical perturbations of the field variables l a p, $, "J on the cosmological spacetime 
manifold in the framework of general relativity. The affine connection on the spacetime manifold of general relativity 
is represented by the Christoffel symbols while the gauge transformation is generated by a flow of an arbitrary vector 
(gauge) field that maps the manifold into itself. Gauge transformation of the fields on a curved manifold is 
associated with a Lie transport of the fields along the vector flow [65, 99]. Infinitesimal gauge transformation is a 
Lie derivative of the field taken at the value of the parameter on the curves of the vector flow equal to 1 [58, chapter 
3.6]. 

Let us consider a mapping of spacetime manifold into itself induced by a vector flow, = £ a (xP). It means that 
each point of the manifold with coordinates x a is mapped to another point with coordinates 



x a = x a -£ a {x) 



(175) 
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This mapping of the manifold into itself can be interpreted as a local diffeomorphism which transforms the field 
variables in accordance to their tensor properties. The transformed value of the field variable is pulled back to the 
point of the manifold having the original coordinates x a , and is compared with the value of the field at this point. 
The difference between the transformed and the original value of the field, generated by the diffeomorphism (175) is 
the gauge transformation of the field that is given by the Lie derivative taken along the vector field at the point 
of the manifold with coordinates x a . 

Let us denote the transformed values of the field variables with a hat. In the linearized perturbation theory of the 
cosmological manifold, the gauge transformations of the field variables are given by equations 

X a (j = X a/j + C a \p + £/3\a ) (176a) 
l a p = la/3 - £a\p ~ €p\a + 9apC h i (176b) 
4> = + *|af , (176C) 

^ = v>+*i«r , (i76d) 

where the hat above each symbol denotes a new value of the field variable after applying the gauge transformation, 
and all functions are calculated at the same value of coordinates x a . The gauge transformations of the field variables 
are expressed in terms of the covariant derivatives on the manifold and, thus, are coordinate-independent. Equation 
(176b) is derived from the Lie transformation (176a) of the metric tensor perturbation, and the relationship (115) 
between x a p and l a p- 

Gauge invariance of the Lagrangian perturbation theory means that the gauge transformations of the field variables 
do not change the content of the theory. In other words, the equations for the field variables must be invariant with 
respect to the gauge transformations (176a)-(176d). However, direct inspection of equations (164), (166), (167) shows 
that they do depend on the choice of the gauge in the form of the gauge function B a introduced in equation (163). 
To find out the gauge-invariant content of the theory one should search for the gauge-invariant field variables and 
to derive the gauge-invariant equations for them. This program has been completed by Bardeen [8] who used the 
functional 3+1 decomposition of the metric tensor perturbations and the vector field £ a to build the gauge-invariant 
variables out of the various projections of the metric tensor components on space an time. Various modifications of 
Bardeen's approach can be found, for example, in [13, 29, 32, 33, 73, 82] and in the book by Mukhanov [81]. We use 
algebraic 3+1 decomposition of the metric tensor perturbations (174) that allows us to build gauge-invariant scalars. 
Vector and tensor perturbations remain gauge-dependent in this approach. In order to suppress the gauge degrees of 
freedom in these variables we impose a particular gauge condition B a — in equation (163). This limits the freedom 
of the gauge field by a particular set of differential equations which are discussed in section (VUG). 



C. The Gauge-invariant Scalars 

The existence of the preferred four- velocity, u a , of the Hubble flow in the expanding universe provides a natural way 
of separating the perturbations of the field variables in scalar, vector, and tensor components. This section discusses 
how to build the gauge-invariant scalars. Vector and tensor perturbations are discussed afterwards. 

The gauge-invariant scalar perturbations can be build from the perturbation of the Clebsch potential, 0, the 
perturbation of the scalar field ip, and two scalars associated with the trace of the metric tensor, q, and its projection 
on the worldlines of the Hubble flow, q. To build the first gauge- invariant scalar, we introduce the scalar perturbations 

Xm = 4- , X q = t- , (177) 

Mm Mq 

that normalize perturbations of the Clebsch potential <j> and that of the scalar field ip to the corresponding background 
values of the specific enthalpy, jl m and fi q . The gauge transformations for the three scalars q, Xm, and x q are obtained 
from (176b)-(176d), and read 

q = q - 2u a u^ aW , (178a) 

Xm = Xm - Uaf , (178b) 
Xq = Xq ~ U a ^ , (178C) 

where we have used the definition of the background four- velocity u a — — <&| a //t m = —^>\ a /p, q in terms of the partial 
derivatives of the background values of the scalar fields $ and Equations (178b), (178c) immediately reveal that 
the linear combination 



X — Xm Xq ! 



(179) 
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is gauge-invariant, \ — Xi that is the diffeomorphism (175) does not change the value of the scalar variable \- 
Two other gauge-invariant scalars are defined by the following equations, 

V m = u a Xm \ a - \ , (180a) 

V q ee u a Xq]a - | , (180b) 

or, more explicitly, 

V m = ^u a ^ a - § + 3^H Xm , (181a) 

Mm ^ C 

"- = ^'"-I +3 *- + ^' (181b) 

where the last terms in the right side of these equations were obtained by making use of thermodynamic relationships 
(18), the equality p q — /2 q , and the equations of continuity (93) and (101) for the density of the ideal fluid, p m , and 
that of the scalar field, p q , respectively. 

One can easily check that both scalars, V m and Vq remain unchanged after making the infinitesimal coordinate 
transformation (175). Indeed, the gauge transformation of the derivatives 

Xm|a = Xm|a ~ HP al j^ - Up^ \ a , (182a) 
Xq|a = Xq|a ~ HP a ^ - Uptf \ a , (182b) 

where P a f) — g a p + UaU/3 is the operator of projection on the hypersurface being orthogonal to the Hubble flow of 
four- velocity u a . Making the coordinate transformation (175), and substituting the transformations of functions q, 
Xm and Xq to the definitions of V m and V q we find 

Vm = V m , Vq = V q , (183) 

that proves the gauge- invariant property of the scalars V m and V q . 

Physical meaning of the gauge-invariant quantity V m can be understood as follows. We consider the perturbation 
of the specific enthalpy p m defined in equation (21). Substituting the decomposition (106) of the field variables to 
equation (21) and expanding, we obtain 

Mm = Am + 5p m , (184) 

where the perturbation Sp m of the specific enthalpy is defined (in the linearized order) by 

Sp m = u a <p\ a - i/2 m q . (185) 



It helps us to recognize that 



V m = ^ + 3^H Xm • (186) 

Mm c 



Fractional perturbation of the specific enthalpy can be re-written with the help of thermodynamic equations (18) in 
terms of the perturbation Se m of the energy density of the ideal fluid, 

^ = 4^!^, (187) 

Mm C ^m "T Pm 

or, by making use of equation (15), in terms of the perturbation 6p m of the density of the ideal fluid 

5p m c 2 Sp m 



Mm C 2 p n 



(188) 



This allows us to write down equation (186) as follows 



V m = £(^+3H Xm ) . (189) 



32 



which elucidates the relationship between the gauge-invariant variable V m and the perturbation Sp m of the rest mass 
density of the ideal fluid. More specifically, V m is an algebraic sum of two scalar functions, Sp m and Xm neither of 
each is gauge- invariant. The gauge transformation of the ideal-fluid density perturbation is 

Spm = Sp m - p m \ a C = Sp m + 3Hp m u a £ a , (190) 

and the gauge transformation of the variable Xm is given by (178b). Their algebraic sum in equation (189) does 
not change under the diffeomorphism (175) showing that V m is the gauge-invariant density fluctuation that does not 
depend on a particular choice of coordinates on spacetime manifold. 

Similar considerations, applied to function V q reveals that it can be represented as an algebraic sum of the pertur- 
bation, Sp q , of the density of the scalar field, and the function x q , 

V q = ^ + 3H Xq . (191) 

It is easy to check that each term in the right side of this equation is not gauge-invariant but their linear combination 
does. The reader should notice that standard textbooks on cosmological theory (see, for example, [71, 84, 99, 100]) 
derive equations for the density perturbations Sp/p but those equations are not gauge-invariant and, hence, their 
solutions should be interpreted with care (see discussion in section IX A 2). 



D. Field Equations for the Scalar Perturbations. 



1. Equation for a scalar q. 



Function q was defined in (168a). In order to derive a differential equation for q, we apply the covariant Laplace- 
Beltrami operator to q, and make use of the covariant equations (161) and (163). Straightforward but fairly long 
calculation yields 



2k 

qi%-2 (if + ff 2 - — I q + 8np m p, 



a- 

~l6np q 



1 K } Vr 



1 + 3^ E X 



(192) 



fdW 

i^ w +2Hfl q )x q ~2u a ^B a ^-4Hu a B a = 8w(a + r) , 



where the source density a + t for the field q is 

o + r = (u a u fj + P afi ) % afj 



(193) 



in accordance with the definitions introduced in (104a), (104b). The reader should notice that equation (192) depends 
on the gauge function B a which remains arbitrary so far. 



2. Equation for a scalar p. 



Function p was defined in (168b). In order to derive equation for p, we apply the covariant Laplace-Beltrami 
operator to the definition of p, and make use of the covariant equations (161) and (163). It results in a wave equation 

4fc 

P'% + ^2 P + B % - 2u a u f3 B alf! - 6Hu a B a = 16ttt . (194) 
where the source density r has been defined in (104b). Equation (194) depends on the arbitrary gauge function B a . 



3. Equation for a scalar x- 



Equation for the gauge invariant scalar, \ — Xm — Xq, is derived from the definitions (177) and the field equations 
(166), (167). Replacing <p an d ?A m those equations with Xm and Xq> an d making use of equations (87), (88) for 
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reshuffling some terms, yields 

xt\ a + 2Hu a Xm \ a "(^-^)xm + W m + (l - ^) u a V m{a - 16^ q X = u a B a , (195a) 
Xq % + 2ffi Q x q | Q - [H - Xq + 4# q + ^fl^q + Wm/imX = ^B a . (195b) 



x'% + 6ffi Q X|« + 3i/x = — ^-Vq - ( 1 - ^ J « a Kn|a • (196) 



Subtracting (195b) from (195a) cancels the gauge-dependent term, u a B a , and brings about the field equation for x, 

Mq9* q I c s 2 

This equation is apparently gauge-invariant since any dependence on the arbitrary gauge function B a disappeared. 
It is also covariant that is valid in any coordinates. 

Equation (196) can be recast to the form of an inhomogeneous wave equation: 

(PmX)'" |a = 2 ^% V * i 1 |) P^ a V m{a . (197) 

Yet another form of equation (196) is obtained in terms of the variable p q x- By simple inspection we can check that 
equation (196) is transformed to 

u d 2 W dW ( c 2 \ 

^ I" - -9*2 M = 2 W Kl ~V~^J P ^ V ^ a ■ (198) 

This is an inhomogeneous Klein-Gordon equation for the field (p q x) governed by V m . The 'mass' of the scalar field 
p q X depends on the second derivative of the potential function W which defines the 'coefficient of elasticity' of the 
background scalar field 

Inhomogeneous equations (196), (197), (198) have the source terms that is determined by variables Vm and V q . We 
derive differential equations for these field variables in the next sections. 



4- Equation for a scalar V m 



Equation for the field variable V m is derived from the equations for functions Xm and q that enter its definition 
(180a). By applying the Laplace-Beltrami operator to function V m we get 



i |c« 



q'% + uPR a pXm\ a + 2Hu a ( V m + -q ) + 3H z [V m + -q 



1 



l» 



1 



(199) 



The Laplace-Beltrami operator for function Xm is given in equation (195a) which is not gauge- invariant. Taking the 
covariant derivative from this equation and contracting it with u a brings about the first term in the right side of 
equation (199), 

„2 ^ 



Ak 



(200) 
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u a X \ a -Gx+ 7 H 1-4 W 



l a u p B. 



a\f3 



The Laplace-Beltrami operator for function q has been derived in (192). Now, we make use of equations (192), (195a), 
(200) in calculating the right side of (199). After a significant amount of algebra, we find out that all terms explicitly 
depending on q and the gauge functions B a cancel out, so that equation for V m becomes 



Vt\ a + (1 



(l - ^ u a u?V mlaP + 2 (3 - J) Hu a V m]a + 



(201) 
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Second-order covariant derivatives in this equation read 



u u r 



|a/3 



:M a # + P 



a/3 



m|a/3 ? 



(202) 



and they form a hyperbolic-type operator describing propagation of sound waves in the expanding universe from the 
source of the sound waves towards the field point with the constant velocity c 2 . Additional terms in the left side of 
equation (201) depend on the Hubble parameter H, and take into account the expansion of the universe. Equation 
(201) contains only gauge-invariant scalars, V m and \- Moreover, it does not depend on the choice of coordinates on 
the background manifold. It also becomes clear that the field variables V m and \ are coupled through the differential 
equations (198) and (201) which should be solved simultaneously in order to determine these variables. Solution of 
the coupled system of differential equations is a very complicated task which cannot be rendered analytically in the 
most general case. Only in some simple cases, the equations can be decoupled. We discuss such cases in section IX. 



5. Equation for a scalar V q . 



The field variable V q is not independent since it relates to V m and x by a simple relationship 



V„ = V m 



X\a 



(203) 



which is obtained after subtraction of equation (180a) from (180b). Equation for V q is derived directly from (203) 
and equations (201) and (196) for V m and x respectively. We obtain, 



W^) U Vqla 



(204) 



21H + 3H 2 
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This equation can be also derived by the procedure being similar to that used in the previous subsection in deriving 
equation for V m . We followed this procedure and confirm that it leads to (204) as expected. Equation (204) is clearly 
gauge- invariant. It couples with the variable \ ancl should be solved along with equation (196). 



E. Field Equations for Vector Perturbations 

Vector perturbations of the ideal fluid and scalar field are gradients, (j)\ a and ip\ a . However, they are insufficient to 
build a gauge- invariant vector perturbation out of the vector perturbation of the metric tensor p a . Field equations 
for vector p a can be derived by applying the covariant Laplace-Beltrami operator to both sides of definition (170) and 
making use of equation (164). After performing the covariant differentiation and a significant amount of algebra, we 
derive the field equation 

pJV - 2Hu a p^ - (^2H + 3H 2 - ^ p a + Pjif (B f3h + B lW + 2Eu 1 B & ) = 16^t q , (205) 

where the matter current a a is defined in (104c). This equation is apparently gauge-dependent as shown by the 
appearance of the gauge function B a . This equation reduces to a much simpler form 

pj^p - 2Hu a pp\f } - (2H + 3H 2 -^p a = 16nr a , (206) 

in a special gauge B a =0 which imposes a restriction on the divergence of the metric tensor perturbation in equation 
(163). 

Equation (205) points out that the vector perturbations are generated by the current of matter r a existing in the 
localized astronomical system which physical origin may be a relict of the primordial perturbations. We do not discuss 
this interesting scenario in the present paper as it would require a non-conservation of entropy and non-isentropic 
background fluid - the case which we have intentionally excluded in order to focus on derivation of cosmological 
generalization of the post-Newtonian equations of relativistic celestial mechanics [58] . 
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F. Field Equations for Tensor Perturbations 

Field equations for traceless tensor p^ can be derived by applying the covariant Laplace-Beltrami operator to the 
definition (171) and making use of equation (164) along with a tedious algebraic transformations. This yields the 
following equation 

P^ l7 |7 - 2ff (u aP ^ + upp^-y) - 2 (tf 2 + A) Pl P^Pf (B, w + B„|„) + \p afi P^B Av = 16^ . (207) 
Here the transverse and traceless tensor source of the tensor perturbations is 

rip = r af 3 - \p a p r , (208) 

where r a p has been introduced in (104d), and r = P afi T a f} in accordance with equation (104b). Tensor is traceless, 
that is g aP r T a0 = P^t^ = 0. 

Equation (207) is gauge-dependent. The gauge freedom is significantly reduced by imposing the gauge condition 
B a = which brings equation (207) to the following form, 



P^ l7 |7 - 2ff(o a p£ 7 lT + u pl^) -2[H 2 + - J pl fj = Wttt^ . (209) 

G. The Residual Gauge Freedom 

The gauge freedom of the theory under discussion is associated with the gauge function B a appearing in equation 
(163). The most favourable choice of the gauge condition is 

B a = , (210) 

which drastically simplifies the above equations for vector and tensor gravitational perturbations. The gauge (210) is a 
generalization of the harmonic (de Donder) gauge condition used in the gravitational wave astronomy and in the post- 
Newtonian dynamics of extended bodies. This choice of the gauge establishes differential relationships between the 
algebraically-independent metric tensor components introduced in section VII A. Indeed, substituting the algebraic 
decomposition (174) of the metric tensor perturbations to equation (163) and imposing the condition (210) yields 

Pali + u a p/3 lf3 + uppj 13 - [u a Ufi - ^Paf^j P l/3 + 2u a u /3 ql /3 + 2Hp a + 2Hqu a = 16n (p m p, m Xm + PqM q x q ) u a . (211) 

Projecting this relationship on the direction of the background 4- velocity, u a , and on the hypersurface being orthogonal 
to it, we derive two algebraically- independent equations between the perturbations of metric tensor components and 
of the matter variables. They are 

p /3 l /3 +%(2q-p) l/3 + 2#q = 16tt (p mJ u m Xm + PqMqXq) , (212a) 
PlpW +* f lpJ f) + lp a ppW + 2Bp a = 0. (212b) 

The gauge (210) does not fix the gauge function £ Q uniquely. The residual gauge freedom is described by the 
gauge transformations that preserve equations (212a), (212b). Substituting the gauge transformation (176b) of the 
gravitational field perturbation l a p to equation (163) and holding on the gauge condition (210), yields the differential 
equation for the vector function £ Q 

r'V + r 1 (^i 7/ 3 - £V) + 2# + i^ a u fi - e^) - ^ (p m ^ m + p^) em a = o , (213) 

which can be further recast to 

+ 2# + t' M u t 3 - S?\pu a ) +2{h-^j fu u a + (h + 3H 2 + C = . (214) 
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The gauge function can be decomposed in time-like, £ = —^up, and space-like, ( a = P a , components, 

Z a = ( a + u a £,. (215) 
Calculating covariant derivatives from £ and C, a and making use of equation (214), yield equations 

^\p + 2Hu%- (h-^)£ -0, (216a) 

C°% + 2H(u?Cw ~ « Q CV) + (h + H 2 + J) C = ' ( 216b ) 

These equations have non-trivial solutions which describe the residual gauge freedom in choosing the coordinates on 
the background manifold subject to the gauge condition (210). It is remarkable that equations (216a), (216b) are 
decoupled and can be solved separately. It means that the residual gauge transformations along the worldlines of 
the Hubble flow are functionally independent of those performed on the hypersurface being orthogonal to the Hubble 
flow. 



VIII. POST-NEWTONIAN FIELD EQUATIONS IN A SPATIALLY-FLAT UNIVERSE 
A. Cosmological Parameters and Scalar Field Potential 

Equations of the field perturbations given in the previous section are generic and valid for any model of the FLRW 
universe. They neither specify the equation of state of matter, nor that of the scalar field, nor the parameter of the 
space curvature k. By choosing a specific model of matter and picking up a value of k = —1,0, +1, we can solve, at 
least, in principle the field equations governing the time evolution of the background cosmological manifold. Realistic 
models of the cosmological matter are rather sophisticated and, as a rule, include several components. It leads to the 
system of the coupled field equations which can be solved only numerically. However, the large scale structure of the 
universe is formed at rather late stages of the cosmological evolution being fairly close to the present epoch. Therefore, 
the study of the impact of cosmological expansion on the post-Newtonian dynamics of isolated astronomical systems 
is based on the recent and present states of the universe. 

Radiometric observations of the relic CMB radiation and photometry of type la supernova explosions reveal that at 
the present epoch the space curvature of the universe, k = 0, and the evolution of the universe is primarily governed 
by the dark energy and dark matter, which make up to 74% and 24% of the total energy density of the universe 
respectively, while 4% of the energy density of the universe belongs to visible matter (baryons), and a tiny fraction of 
the energy density occupies by the CMBR radiation [35, 47, 53, 55]. It means that we can neglect the effects of the 
baryonic matter and CMB radiation field in consideration of the post-Newtonian dynamics of astronomical systems 
in the expanding universe. 

We shall assume that the dark matter is made of an ideal fluid and the dark energy is represented by a scalar field 
with a potential function W which structure should be further specified. In doing this, we shall follow discussion in 
[2] assuming that the spatial curvature k = 0, and the potential, W, of the scalar field relates to its derivative by a 
simple equation 

dW , — - 

^ = -y/to\W , (217) 

where the time-dependent parameter, A = A(^>), characterizes the slope of the field potential W. The time evolution 
of the background universe can be described in terms of the parameter A and two other parameters, xi = xi(<I>) and 
x 2 = x 2(5 r ): which are functions of the density, p q = j2 q = ^, of the background scalar field, and the potential, W, 
scaled to the Hubble parameter, H. These parameters are defined more specifically as follows, 

P1 = ^ X1 , (218) 

W = ^-x 2 . (219) 

87T 

The energy density of the scalar field, e q , is expressed in terms of the parameters xi and x 2 and the parameter 
fi q = 87re q /3-ff 2 , by a simple relationship 

Ci q = xi + x 2 . (220) 
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The time evolution of the parameters xi and x 2 is given by the system of two ordinary differential equations which 
are obtained by differentiating the definitions (218), (219) and making use of the equations (101) along with the 
Friedmann equation (87) with k = 0. It yields 

— - = -6x 1 + AV6x7x 2 + 3x 1 [(l-w m )xi + (l + w; m )(l-x 2 )] , (221a) 

^ = -AV6x7x 2 + 3x 2 [(1 - w m ) Xl + (1 + w m ) (1 - x 2 )] , (221b) 

where uj = In a is the logarithmic scale factor characterizing the number of e-folding of the universe, w m is the 
parameter entering the hydrodynamic equation of state (89), and the parameters xi and x 2 are restricted by the 
condition imposed by the Friedmann equation (85), that is 



xi+x 2 = l-fi m , (222) 



where f2 m = 87re m /3-ff 2 . 

The parameter A obeys the following equation 

dX 

duj 

where 



6 Xl A^ (r q - 1) , (223) 



T q = '- o W , 224 

q (dw/d*) 2 

If T q = 1, the parameter A is constant, and equation (217) can be integrated yielding an exponential potential 

W(V) = W exp(- V8^A§) . (225) 

In this case, and under assumption that, w m — const., the system of two differential equations (221a), (221b) is 
closed. If T q 7^ 1, three equations (221a)-(223) must be solved in order to describe the evolution of the background 
cosmological manifold. 

In the general case, derivatives of the potential W are expressed in terms of the parameters under discussion. 
Namely, 



9W _ 3A 2 d 2 W _ x2o2 



ff 2 x 2 , _ = SrqA 2 ^ . (226) 



It is also useful to express the products p q ^ q and p m ^ m in terms of the parameters xi and x 2 . For ^ q = p q , one can 
use definition (218) to obtain 

3H 2 

h^q. = xi . (227) 

The product p m ji m = e m +p m , so that making use of the matter equation of state, p m = w m e m , and equation (222), 
we derive 

3H 2 

PmMm = -5— (1 + W m )n m , (228) 

where f2 m = 1 — xi — x 2 . These equations allow us to recast equation (87) for the time derivative of the Hubble 
parameter to the following form 

H = ~(l + w eS )H 2 , (229) 

where 

W c ff = Wm + (1 - tu m )xi - (1 + w m )x 2 , (230) 

is the (time-dependent) parameter of the effective equation of state of the mixture of the ideal fluid and the scalar 
field. 
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B. Conformal Cosmological Perturbations 



The FLRW metric (52) is a product of the scale factor a and a conformal metric f a p- The conformal spacetime 
is comoving with the Hubble flow and is not globally expanding. In case of the flat spatial curvature, k — 0, the 
conformal spacetime becomes equivalent to the Minkowski spacetime which is used as a starting point in the standard 
theory of the post-Newtonian or post-Minkowskian approximations [23]. Therefore, it is instructive to formulate the 
field equations for cosmological perturbations in the conformal spacetime. 

Let us define the cosmological perturbations, /i Q( g, of gravitational field in the conformal spacetime with the back- 
ground metric f a p as follows, 

x a fi = a 2 (v) h ap , (231) 

where perturbations x a p has been defined in (106), and a(r]) is the scale factor of the FLRW universe. Perturbation 
l a p relates to x a p by equation (115), and can be also represented in the conformal form 

Lp = a 2 (r])j a p , (232) 

where 

lap = -/la/3 + aph , (233) 

with h = f a ^h a p. In what follows, tensor indices of geometric objects in the conformal spacetime are raised and 
lowered with the help of the conformal metric f a p. 

We assume that the scale factor a of the universe remains unperturbed. This assumption is justified since we can 
always include the perturbation of the scale factor to that of the conformal metric. Thus, the perturbed physical 
spacetime interval, ds, of the FLRW universe relates to the perturbed conformal spacetime interval, ds, by the 
conformal transformation 

ds 2 = a 2 {ff)ds 2 . (234) 

The perturbed conformal spacetime interval reads 

ds 2 = f a pdx a dx' 3 , (235) 

where 



.fa/3 = f a/3 + h a f) , (236) 

is the perturbed conformal metric. Here, f a p is the unperturbed conformal metric defined in (53), h a p is the per- 
turbation of the conformal metric, and x a — (x°,x l ) are arbitrary coordinates which are the same as in the physical 
spacetime manifold. 

It is worth emphasizing that in case of the space curvature k — 0, the background conformal metric, g a p(r],X l ), 
expressed in the isotropic Cartesian coordinates (ij, X 1 ). is the diagonal Minkowski metric, g Q , ( g(?7,X t ) = rj a p = 
diag(— 1, 1, 1, 1). Therefore, the background metric f a p remains the Minkowski metric with the components expressed 
in arbitrary coordinates by means of tensor transformation 

f a p = M^M"^ - ( 23? ) 

where the matrix of transformation has been defined in (50). If the matrix of transformation, M.^ a , is the Lorentz 
boost, the conformal metric, f a p, remains flat, f a p = rj a p. It is worth noticing that, in general, the unperturbed 
conformal metric can be chosen flat even in case of k = — 1,+1 [49]. It means that our formalism is applicable to 
FLRW universe with any space curvature. However, the conformal factor in this case is not merely the scale factor 
a(ri) of the FLRW universe but a more complicated function of coordinates. Though it is not difficult to handle all 
three cases of k = — 1, 0, +1 but it burdens equations for the field perturbations and we restrict ourselves only to the 
case of the spatially flat universe with k = 0. 

Similarly to (174) the conformal metric perturbation, j a p 7 can be split in the algebraically-irreducible components 

lap = P T a p + V a pp + V/3p a + (v a v p + l-n a p J p + 2v Q v ) g (q - p) , (238) 
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where the four-velocity v" = au a , v a = f a p^ — a 1 g a /3U^ — a 1 u a . and 

= f a/3 + V a V/3 , (239) 

is the operator of projection on the conformal space which represents a hypersurface being everywhere orthogonal to 
the congruence of worldlines of four- velocity v a . Four- velocity v Q is an analogue of the Hubble flow in the conformal 
spacetime. We also notice that P a p = a 2 Tt a p. 

Different pieces of the conformal metric perturbation, 7 Q( g, are related to those of the physical metric perturbation, 
la/3, by the powers of the scale factor, 



More specifically, 



Pip = a2 P J a p » Pa = ap a , P=P, q = q ■ (240) 



q = I(W+^)TV, (241a) 

P = if'V , (241b) 

Pa = -7f Q ,3 v 7 7 ( 9 7 , (241c) 

Paf) = Pa/3 ~ ^C0P , (241d) 



where 



Pa/3 = K^TT^l^ ■ (242) 

The trace of the gravitational perturbation, 7 = f a ^'j a p = 2(p — q). The components h a p = —~f a [j + f q^7/2 are used 
in calculating dynamical behavior of particles and light in the conformal spacetime as well as in matching theory with 
observables. The components of h a p are 

2_ 

h a p = -P 7 a p - V a p f3 - V f:S p a + -TT a fiP - (v a Vf3 + TT afj ) q , (243) 

and h = i af3 h af s = 2(p - q) = 7. 

It turns out that the conformal Hubble parameter, !K = a' / a is more convenient in the conformal spacetime than 
H = R/R = R~ 1 dR/dT, where T is the cosmological time (see section VB). Relationships between J£ and H, and their 
derivatives are shown in equations (42)-(44). These relationships along with equations (43) and (229) are employed 
in order to express the time derivative, 5f', of the conformal Hubble parameter in terms of Jf 2 and the parameter w e s 
of the effective equation of state 

W = -i(l + 3 Wcff )M 2 . (244) 
We shall use this expression in the calculations that follows. 



C. The Post-Newtonian Field Equations in Conformal Spacetime 



The set of the post- Newtonian field equations in cosmology consists of equations for perturbations of the background 
matter and gravitational field. Perturbations of matter are described by four scalars, V m , V q , \m and x q but only 
three of them are functionally- independent because of equality (203), that is 

V m - l/ q = U a ( Xm - Xq)| Q ■ (245) 

Depending on a particular situation, any of the three scalars can be taken as independent variables. 

The gravitational field perturbations are q, p, p a , p T aj3 but among them the scalar q is not independent and can be 
expressed either in terms of Xm and V m in accordance with (180a), 

q = -2(V m -u a Xm<a ) , (246) 

where we have also used the equality q = q as follows from (240). The scalar q can be also expressed in terms of 
X q and V q in accordance with (180b). Hence, as soon as the pairs, V m and Xm or V q and x q are known, the scalar 
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gravitational perturbation q can be easily calculated from (246). Functions p, p a , p T a/3 are independent and decouple 
both from each other and from the other perturbations. Thus, the most difficult part of the theory is to find out 
solutions of the scalar perturbations which are coupled one to another. 

The field equations in the conformal spacetime for variables Xm, Xq, V m and for p, p a , p J af3 are derived from the 
equations of the previous section by transforming all functions and operators from physical to conformal spacetime. 
The important part of the transformation technique is based on formulas converting the covariant Laplace-Beltrami 
wave operators, defined on the background spacetime manifold, to their conformal spacetime counterparts. 



1. The Laplace-Beltrami operator 



Let F be an arbitrary scalar, F a - an arbitrary covector, and F a p - an arbitrary covariant tensor of the second 
rank. We have three Laplace-Beltrami operators on the curved background manifold: scalar - F^^, vector - i r c J'*| M , 
and tensor - F a p^\^ types where the covariant derivatives are taken with respect to the affine connection f Q ^ 7 being 
compatible with the metric g a p (see equation 58). Covariant derivatives are the most convenient for the invariant 
description of differential equations of mathematical physics on curved manifolds. For practical purposes of finding 
solutions of the differential equations, the covariant operators must be expressed in terms of partial derivatives with 
respect to coordinates chosen for solving the equations. 

Transformation of the covariant Laplace-Beltrami operators to the partial derivatives is achieved after writing down 
the covariant derivatives for scalar, vector and tensor in explicit form by making use of the Christoffel symbols given 
in (60)-(62). Tedious but straightforward calculations of the covariant derivatives yield the scalar, vector and tensor 
Laplace-Beltrami operators in the following form 



!" a 2 



DF - 2 < Kv< x F. ll 



F K 



1 



F aK — 



l r 

a 



DF a - 2Xv"F ll . a + 2'Kv a t liv F jl . v + [pd + 2'K 2 ) F a - 2X 2 v a v"F tl 



(247a) 
(247b) 

DF afj + 2Mv»F a ^ - m^F^ - 25Cv"F„ j8 ; a + 2X1 » v {v a F^, v + Y fi F a ^, v ) (247c) 
+ 2 (5C' + X 2 ) F a /j - m 2 (>v a F^ + v»Vf,F ail - U^^F^ - if a ^v v F^~ 
where we have introduced notations 



DF = f^F, 



DF„ 



OF, 



a/3 



(248) 



of the wave operators for scalar, vector and tensor fields in the conformal spacetime and in arbitrary coordinates. 
Notice that although the conformal spacetime coincides, in case of k = 0, with the Minkowski spacetime, the metric 
f a/3 is not the diagonal Minkowski metric r) a p unless the coordinates are Cartesian. Of course, the covariant derivative 
from a scalar must be understood as a partial derivative, that is F. a — F„. 

We will need several other equations to complete the transformation of the Laplace-Beltrami operators to the 
conformal spacetime since the wave operator □ acts on functions like (240) which are made of a product of some 
power n of the scale factor, a — 0(77), with a geometric object, F = F{x a ), which can be a scalar, a vector or a tensor 
of the second rank (we have suppressed the tensor indices of F since they do not interfere with the derivation of the 
equations which follow). These equations are 

(a n F).„ = a" (F;„ - nXv^F) , (249a) 
(aY).^ = a" [F ;M „ - WK (v M F ;v + v„F ; „) + n (5£' + nX 2 ) v„v„] , (249b) 

and they allow us to write down the wave operator from the product of a" and F in the following form 



□ ( a n F) = a n IDF - 2nJ£v Ai F ;AJ - n (?f + WK 2 ) F 



(250) 



It is easy to confirm that contraction of (249b) with the conformal four- velocity, v Q , brings about another differential 
operator 



(a n F).^ = a n ^v v F;„ v + 2nXv»F^ + n (JC' + WK 2 ) F 



(251) 
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We remind that if the object f is a scalar, the covariant derivative is reduced to a partial derivative, F - a = F , a - In 
case, when F is either a vector or a tensor, the covariant derivative must be calculated with taking into account the 
affine connection B a ( 9 7 defined in (62). 

It is also interesting to notice that in the expanding universe the conformal Laplace operator, AF = n^F-^v is the 
scale invariant in the sense that 

A (a n F) = a n AF , (252) 

where F is a tensor of an arbitrary rank. Equation (252) can be proven by adding up (250) and (251), and accounting 
for definition (239) of the projection operator on the hypersurface being orthogonal to v°. 

Now, we are ready to formulate the field equations for cosmological perturbations in the conformal spacetime. 



2. Equations for the matter perturbations 

We accept the gauge condition imposed by equations (163), (210) and convert the covariant derivatives taken with 
respect to the background metric, g a /3, to the partial derivatives of the conformally-flat metric, f a pi m equation (201) 
for scalar V m . We use equation (247a) for the Laplace-Beltrami operator, and expressions for various cosmological 
parameters given in section VIII B. After arranging terms with respect to the powers of the Hubble parameter J£, we 
obtain the sound-wave equation for function V m describing perturbations of the ideal fluid, 



(253) 



s 1 - w cB - i(i + Wm ) n - ^ 1 <>„ 



+ 12Jf 2 



ft v m 




— = -Ana 2 (a + r) 



Similar procedure applied to equation (204) leads to a wave equation for function V q describing perturbations of the 
scalar field, 



DV q + 2[l- x /— Ax 2 ) Mv^ q , M 



(254) 



+3 



1 - w cH - -(1 + w m ) ( 1 



+Ax 2 



3A 2r, 



+hi 2 {i + Wm ) (3 



X2 
XI 



v"x,m - 3-| 



x 2 v q 

a 



= -Ana 2 (a + r) 



Equations (253) and (254) contains function x which obeys equation (196). Making use of the same transformations 
as above, we recast (196) to a wave equation for \, 



□X + 4x(l- y^^) V«X,« - I (1 + Wef) = -O 



— AxX.Yl , • ( 1 '— ) v'M;,,,, 



(255) 



Equations (253)-(255) are closed with respect to the variables V m , V q and x- The gauge-invariant scalar, V m describes 
propagation of sound waves in the ideal fluid filling up the expanding universe. It can be found from solving two 
equations (253) and (255) simultaneously after imposing a certain (cosmological) boundary conditions. As soon as the 
gauge- invariant scalar \ is known, the potential, V q , can be determined as a particular solution of the inhomogeneous 
equation (254) or, more simple, from equation (203). 

We also need equations for the normalized Clebsch and scalar potentials, Xm and Xq- These potentials are required 
to determine the gravitational perturbation, q, with the help of (246) and/or to get the check on self-consistency of 
the solutions of equations in the matter sector of the perturbation theory. Conformal-spacetime equations for Xm and 
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X q are derived from their definition (177) and the field equations (166) and (167). They are 



□Xm + 2 (1 + w cff) ^ 2 Xm = 121K 2 x lX - a ^KV™ + \l - - J v a V m . 
□Xq + | (1 + Weff) W 2 Xq = -6W 2 (1 + Wm )fi m x -a (a- "KVq . 

By subtracting one of these equations from another, we get back to equation (255). 



(256) 
(257) 



3. Equations for the metric perturbations 

Post-Newtonian equations for gravitational perturbations in physical spacetime are (192), (194), (205) and (207). 
We remind to the reader that the gauge conditions (163), (210) has been imposed. In this gauge, equations for the 
conformal metric tensor perturbations become 



Dq - 2Kv a q ia + (1 + 3w eS ) = 87ra 2 (a + r) - 24Jf 2 



-Ax2 — "Kx\ 



Xq 

a 



Up - 21Kv> jQ 
D Pa - 25tv /3 p Q;/3 + (1 + 3w cfi ) l K 2 p a 
Uplp - 25CvTp aj8 ; 7 



-3 (1 + w cS ) K 2 n n 

167ra 2 r , 
16tt ar a , 

167TT 



1-3)^-^1 + 3^ 



2 \ 
c = \ Xn 



afj ■ 



(258a) 



(258b) 
(258c) 
(258d) 



The reader can observe that equations (258a)-(258d) for linearized metric perturbations are decoupled from each 
other. Moreover, equations (258b-(258d) are decoupled from the matter perturbations V m , Xm, etc. Only equation 
(258a) for q is coupled with the matter perturbations governed by equations (253), (256), (257) so that these equations 
should be solved together. As we have mentioned above, function g is a linear combination of V m and Xm according to 
(246). Hence, in order to determine q it is, in fact, sufficient to solve (253) and (256). Nevertheless, it is convenient to 
present the differential equation (258a) for q explicitly for the sake of mathematical completeness and rigour. It can be 
used for independent validation of the solution of the system of equations (253), (256) and (246). Unfortunately, these 
equations are strongly mixed up and cannot be solved analytically in the most general situation of a multi-component 
background universe governed by the dark energy and dark matter. Solution of (253)-(257) and (258a)-(258d) would 
require an application of the methods of numerical integration. 

It would be instrumental to get better insight to the post-Newtonian theory of cosmological perturbations by making 
some simplifying assumptions about the background model of the expanding universe in order to decouple the system 
of the post-Newtonian equations and to find their analytic solution explicitly. We discuss these assumptions and the 
corresponding system of the decoupled post-Newtonian equations in the section IX. 



D. The Residual Gauge Freedom in the Conformal Spacetime 

The gauge conditions (163), (210) in the physical space are given by (212a), (212b). After transforming to the 
conformal spacetime the gauge conditions read 

P 13 ;P + v fj (2q - p) fj + 23{q = 16™ (p m p m Xm + PqMqXq) . (259a) 

p ra ^. i p+v p p a .,p + \w ap p,p + 2Kp a = 0. (259b) 

The residual gauge freedom in the conformal spacetime is described by two functions, £ = £/a and ( a . where £ and 
( a have been defined in section VUG. Differential equations for £ and ( a are obtained by making transformation of 
equations (216a), (216b) to the conformal spacetime. The calculation is straightforward and results in 



□C - 2!Kv /3 C/3 + (1 + 3w cS ) J£ 2 C = 0, 
□C a -2Mv /3 C Q ; ^ = 0. 



(260a) 
(260b) 
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Solutions of equations (258a)-(258d) are determined up to the gauge transformations 

q = g + 2v a C,a + 2MC, (261a) 

P = P + C% + 3v Q Ca + 65£C , (261b) 

p a = p Q +^(vTC ,3 ;7-C" 3 + 2^) , (261c) 

p af 3 = P a fi ~ (TT^aV + Tf^Tfa") CV + * af 3 (C ;a + V^C.a + 20£C) , (261d) 

where the gauge functions (, ( a are solutions of the differential equations (260a), (260b). 



IX. THE DECOUPLED SYSTEMS OF THE POST-NEWTONIAN FIELD EQUATIONS 
A. The Universe Governed by the Ideal Fluid and Cosmological Constant 

1. Case 1: Arbitrary equation of state of the ideal fluid 

Let us consider a special case of the dark energy represented by the cosmological constant A. In this case, the 
equation of state of the scalar field is w q — 1, and we have p q p q — e q + p q = 0. The parameter xi = 0, and 
x 2 = A/(3-ff 2 ). It yields the parameter f2 q = x 2 , and f2 m = 1 — x 2 . Since the cosmological constant corresponds to a 
constant potential W of the scalar field, we get for its derivative dW /d^> = 0, and equation (217) points out that the 
parameter A = 0. 

In the universe governed by the ideal fluid and the cosmological constant the parameter of the effective equation of 
state 



W e ff = w m - (1 + fflm)-i • (262) 

Oil 



Hence, the time derivative of the Hubble parameter defined in (229), is reduced to a more simple expression, 

H= i(l + ™ m )(A-3ff 2 ) . (263) 
On the other hand, equation (87) tells us that in this model of the universe 

H = -4irp m fL m , (264) 

The field equation (253) for scalar V m is reduced to that describing the time evolution of the perturbation of the ideal 
fluid density, 8p m . Indeed, the scalar V m defined by equation (180a), can be recast to the form given by equation 
(189), that is 

V m = ^Sm, (265) 



where the gauge-invariant scalar perturbation 



S m = + oH Xm , (266) 

Pm 

is a linear combination of the perturbation of the mass density of the fluid and the normalized Clebsch potential. 
Replacing expression (265) in equation (253), yields the exact equation for 6 m that is 

2 \ 2 / 2 \ 

v a v^ m;a/3 - C -^U8 m + h - 3^J Xv a 8 m , a (267) 



c 2 

(1 - 3w m ) -§ + (1 + W m ) 



K Sm + (1 + w m ) l-3-| a'AS m = 4ira 2 (a + r) . 



This equation describes propagation of the ideal fluid density perturbation 6 m in the form of sound waves with velocity 
c s - 
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Equation (255) for potential \ makes no sense since the normalized perturbation Xq = i>/fiq of the scalar field 
diverges due to the condition fi q = p q = 0. Equation for the perturbation of the scalar field tp itself is obtained from 
(167) and is reduced to a homogeneous wave equation 

□V> - 25^^ = . (268) 

Equation for the normalized Clebsch potential, Xm, is derived from equation (256). In the case of the universe 
under consideration this equation reads 

□Xm + \ (1 + w m ) (3J£ 2 - a 2 A) X m = (l - av"$ mi/1 - Aa < K < ^5 m . (269) 

This is an inhomogeneous equation that can be solved as soon as one knows <5 m from equation (267). 

Gravitational potential q can be determined directly from equation (246) after solving equations (267) and (269) 
or by solving equation (258a) which takes on the following form, 

,2\ / c 2\ 

c s \ Xn 



Dq - 2Jfv^ + [(1 + 3w m ) IK 2 - (1 + w m ) a 2 A] q = 8na 2 |a + r + p m p m (l - ^ «5 m + "K (l - 3 , 



Cr I a 



(270a) 



Equations for the remaining gravitational perturbations are found from (258b)-(258d) which read 

Dp - 2%v^p^ = 167ra 2 T , (270b) 
Up a - 2^Kv^p a;ft + [(1 + 3w m )^ 2 - (1 + w m )a 2 A] p a = \^ar a , (270c) 

Upl p 25C^. M = 16^ . (270d) 

Equations given in this section are valid for arbitrary cosmological equation of state of the ideal fluid, p m — ui m 6 m , 
that is physically reasonable. The parameter w m of the equation of state should not be replaced with the ratio of 
c 2 /c 2 which characterizes the derivative of pressure p m with respect to the energy density e. This is because the 
parameter w m can depend in the most general case on the other thermodynamic quantities (like enthropy, etc.) which 
may implicitly depend on I. Equations (267)-(270d) are decoupled in the sense that all of them can be solved one 
after another starting from solving equation (267) for 5 m , which is a primary equation. 



2. Case 2: Dust equation of state 



Equations of the previous section can be further simplified for some particular equations of state of the ideal fluid. 
For example, in the case when the ideal fluid is made of dust, the background pressure of matter drops to zero making 
parameter of the equation of state w m = 0. Sound waves do not propagate in dust. Hence, the speed of sound 
c s = 0. For this reason all terms being proportional to c 2 and w m vanish in equation (267). Moreover, dust has 
the specific enthalpy, /i m = 1 making the energy density of dust equal to its rest mass density e m = p m , and the 
normalized perturbation Xm of the Clebsch potential of dust is equal to the perturbation <f> of the Clebsch potential 
itself, Xm — </>■ The Friedmann equation (85) tells us that 

= y (87rp m + A) . (271) 

Accounting for this result in equation (267), and neglecting all terms being proportional to the speed of sound, c s , we 
obtain 

v a v^ m;a|9 + < Kv a 8 m ^ a - 4Tra 2 p m 5 m = 4na 2 (a + r) , (272) 

where the terms depending on the cosmological constant, A, have cancelled out. This equation is more familiar when 
is written down in the preferred FLRW frame, where v Q = (1,0,0,0). Equation (272) assumes the "canonical" form 

S m + "K5 m - 4na 2 p m 5 m = 4-ko 2 (a + r) , (273) 

which can be found in many textbooks on cosmology [71, 81, 84, 99, 100]. 

Equation (273) has been derived by previous researchers without resorting to the concept of the Clebsch potential 
of the ideal fluid. For this reason, the density contrast, <5 m , was interpreted as the ratio of the perturbation of the 
dust density to its background value, S — 5p m /p m , without taking into account the perturbation, <ft, of the Clebsch 
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potential. However, the quantity 5 is not gauge-invariant which was considered as a drawback. The scrutiny analysis 
of the underlying principles of hydrodynamics in the expanding universe given in the present paper, reveals that 
equation (273) is, in fact, valid for the gauge-invariant density perturbation <5 m defined above in (266). Another 
distinctive feature of equation (273) is the presence of the source of a bare perturbation in its right side. The bare 
perturbation is caused by the effective density a + r of the matter which comprises the isolated astronomical system 
and initiates the growth of instability in the cosmological matter that, in its own turn, induces formation of the large 
scale structure of the universe [84, 100]. Standard approach to cosmological perturbation theory always set a + t = 
and operates with the spectrum of the primordial perturbation of the density Sp m / p m (but not with the spectrum for 

Equation (269) in case of dust reads, 

□Xm + \ (3^ 2 - a 2 A) X m = av a S m . a , (274) 

where % m is reduced to the perturbation of the Clebsch potential, % m = <f>, for the reason that has been mentioned 
above. 

If equations (273) and (274) are solved, the gravitational perturbations can be found from equations (270a)-(270d), 
which take on the following form 

Uq - 2'Kv a q. a + (Jf 2 - a 2 A) q = 8ira 2 a + T + p m (s m + , (275a) 

Dp - 2JCv a p, Q = 167ra 2 T , (275b) 
Up a - 2Wv f3 p a:JS + (J{ 2 - a 2 A) p a = 16TTaT a , (275c) 
Uplp 2J£v>^ ;7 = ltorr^ . (275d) 

It is interesting to notice that besides the bare density perturbation, a + r, the source for the scalar gravitational 
perturbation, q, contains the induced density perturbation p m (<5 m + 5fx m /a) = Sp m + Hp m <p m the right side of 
equation (275a). This induced density perturbation changes the initial mass of the isolated astronomical system in 
the course of the evolution (expansion) of the universe. This explains the origin of the time-dependence of the central 
point-like mass in the cosmological solution found by McVittie [78] (see also discussion in [16]). 



B. The Universe Governed by a Scalar Field 



In this section we explore the case of the universe governed primarily by a scalar field with all other matter variables 
being unimportant. In this case, the time evolution of the background universe is defined exceptionally by equations 
(221a), (221b). The most general solution of (221a), (221b) is complicated and can not be achieved analytically. 
Numerical analysis shows that the solution evolves in the phase space of the two variables {xi,X2} from an unstable 
to a stable fixed point by passing through a saddle point [2]. The cosmic acceleration is realized by the stable point 
with the values of xi = A 2 /6 and x 2 = 1 — A 2 /6, which is equivalent to the equations of state (89) with the values of 
the parameters, w m — 0, and, w q = — 1 + A 2 /3. It also requires the energy density of the background matter e m = 0, 
that is fi m = 0. In such a universe the derivatives of the potential of the scalar field are 

Moreover, because p m ^ m = e m + p m = 0, the time derivative of the Hubble parameter is 

H = -4^ q = -^H 2 (1 + w q ) . (277) 

In the point of the attractor of the scalar field, perturbations of the ideal fluid are fully suppressed that is the 
Clebsch potential of the fluid, Xm — 0. It makes the function V m = q/2, that is reduced to the perturbation of the 
scalar component of the gravitational field only. Perturbations of the scalar field are described by the scalar field 
variable, x q . In particular, after substituting the derivatives (276) of the scalar field potential along with the derivative 
(277) of the Hubble parameter, in equation (204), one obtains the post-Newtonian equation for function Vq, 

nV q - (1 - 3w q ) Mv"y q ,, + ^K 2 (1 - w q ) (1 + 3w q ) V q = -47ra 2 (<r + r) . (278) 
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Field equation for the perturbation of the scalar field, % q , is reduced to 

□x q - 2Mv^Xq,M + (1 + 3w q ) Xq = - (1 + 3w q ) KV q , (279) 

where the variable, \q = Xq/ a i has been used for the notational convenience. 

Post- Newtonian equations for gravitational perturbations are (258a)-(258d). After substituting the values of the 
parameters xi,x 2 , w c s, etc., corresponding to the model of the universe governed by the scalar field alone, the post- 



Newtonian equations for the metric perturbations become 

Uq - 2ftv fl q^ + (1 + 3w q ) ft 2 q = 87m 2 (a + r) + 3 (1 + w q ) (1 + 3w q ) V?x q , (280a) 

Up - 25Cv^, M = 167ra 2 T , (280b) 

□p Q -2Jfv> Q;AI + (l + 3w q )5{ 2 p Q = 167raT a , (280c) 

□p^ - mVpl^ = 16^ . (280d) 



One can see that the field equations for the scalar field and metric perturbations are decoupled, and can be solved 
separately starting from the primary equation (278). 



C. Post-Newtonian Potentials in the Linearized Hubble Approximation 



1. The metric tensor perturbations 



The post-Newtonian equations for cosmological perturbations of gravitational and matter field variables crucially 
depend on the equation of state of the matter fields in the background universe. It determines the time evolution of 
the scale factor a = a(rj) and the Hubble parameter "K = Ji(r/) which are described by the wide range of elementary 
and special functions of mathematical physics (see, for example, the books by Amendola and Tsujikawa [2], Macias, 
A., Cervantes-Cota, J. L. & Lammerzahl, C. [75], Stephani et al. [94] and references therein). It is not the goal of the 
present article to provide the reader with an exhaustive list of the mathematical solutions of the perturbed equations 
which requires a meticulous development of cosmological Green's function (see, for example, [45, 68, 69, 86]). 

In this section we shall focus on the observation that the post-Newtonian equations for the field perturbations 
have identical mathematical structure if all terms that are quadratic with respect to the Hubble parameter, "K, are 
neglected. In such a linearized Hubble approximation the differential equations for cosmological perturbations are 
not only decoupled from one another, but their generic solution can be found irrespectively of the equation of state 
governing the background universe. Indeed, if we neglect all quadratic with respect to "K terms, the field equations 
for the conformal metric perturbations are reduced to the following set, 

Uq-2'Kv a q. a = 87ra 2 (cr + r) , (281a) 

Up - 2'Kv a p. a = 167ra 2 r , (281b) 

Up a - 2'Kv p a ^ = l6irciT a , (281c) 

Upl 25£vX ft7 = , (281d) 

where the wave operator □ has been defined in (248), and the source of the perturbation is the tensor of energy- 
momentum of a localized astronomical system with the matter having a bounded support in space - see section (VH). 
The differential structure of the left side of equations (281a)-(281d) is the same for all functions. The equations differ 
from each other only in terms of the order of !K 2 which have been omitted. 

In order to bring equations (281a)-(281d) to a solvable form, we resort to relationship (250) which reveals that in 
the linearized Hubble approximation, the equations can be reduced to the form of a wave equation 

U(aq) = 87ra 3 (ct + t) , (282a) 
U(ap) = 16ira 3 T , (282b) 
□ (ap a ) = 16ira 2 T a , (282c) 

° ( a P T a p) = 16*07-^ . (282d) 

So far, we did not impose any limitations on the curvature of space that can take three values: k = {— 1,0, +1}. 
Solution of wave equations (282a)-(282d) can be given in terms of special functions in case of the Riemann (k = +1) 
or the Lobachevsky (k = —1) geometry [68, 69]. The case of the spatial Euclidean geometry (k = 0) is more 
manageable, and will be discussed below. 
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If the FLRW universe is spatially-flat universe, k = 0, and we chose the Cartesian coordinates x a related to the 
isotropic coordinates X a of the FLRW universe by a Lorentz transformation, the operator □ becomes a wave operator 
in the Minkowski spacetime, 

□ = • (283) 

In this case, equations (282a)- (282d) are reduced to the inhomogeneous wave equations which solution depends 
essentially on the boundary conditions imposed on the metric tensor perturbations at conformal past-null infinity 0~ 
of the cosmological manifold [79]. We shall assume a no- incoming radiation condition also known as Fock-Sommerfeld's 
condition [23, 36] 

lim n 7 <9 7 [a(r/)rZ Q/3 (x 7 )] = , (284) 

t-f r=const. 

where x 1 — (x°,x l ), r\ = ?7(x 7 ), the null vector n a = {l,x l /r}, and r = 8i^x l x^ is the radial distance. This condition 
ensures that there is no infalling gravitational radiation arriving to the localized astronomical system from the future 
null infinity 3 + ■ Effectively, it singles out the retarded solution of the wave equation. 

A particular solution of the wave equations satisfying condition (284), is the retarded integral [65] 

aq(t,x) = -2 [ * 3 lv(s^)]l«(s, X i)+T(s, X i)}d3x< ^ 

Jv \ x x I 

. f a 3 \ti(s,x')]t (s 7 x') d 3 x' 
ap(t,x) = -4 / — uy ' n v — '- , (285b) 



\x — x'\ 




a 2 [r](s,x')] r a (s, 


x') d 3 x' 


\x — x'\ 




a [r](s, x')} (s. 


x') d 3 x' 


\x — x'\ 



ap a (t,x) = -4 / L,v ' , /J °V ' , (285c) 

Jv \x — X \ 

f a \ri(s, x')] t! r (s, x') d 3 x' 
apUt,x) = -4 / . (285d) 



where the scale factor a in the left side of all equations is a = a [rj(s, x)], and the argument s of the functions appearing 
in the integrands, is the retarded time 

s = t - \x - x'\ . (286) 

The retarded time s is a characteristic of the null cone in the conformal Minkowski spacetime that determines the 
causal nature of the gravitational field of the localized astronomical system in the expanding universe with k = [58] . 
Solutions (285a)-(285d) are Lorentz-invariant as shown by calculations in Appendix A. 

Integration in (285a)-(285d) is performed over the volume, V, occupied by the matter of the localized astronomical 
system. In case of the system comprised of N massive bodies that are separated by distances being much larger 
than their characteristic size, the matter occupies the volumes of the bodies. In this case the integration in equations 
(285a)-(285d) is practically performed over the volumes of the bodies. It means that each post-Newtonian potential 
q,p,Pa,P T a i3 is split in the algebraic sum of TV pieces 

W JV N N 

q = ^2<iA, p^^pa, p a = ^p Aa , p«/3 = X]pL/3> ( 287 ) 

A=l A=l A = l A=l 

where each function with sub-index A has the same form as one of the corresponding equations (285a)-(285d) with 
the integration performed over the volume, V A , of the body A. 

2. The gauge functions 

The residual gauge freedom describe the arbitrariness in adding solution of homogeneous equations (285a)-(285d) 
with the right side being equal to zero. It is described by two functions, ( = £/a and ( a . Since we neglected the 
terms being quadratic with respect to the Hubble parameter, the gauge functions satisfy the following equations 

□C - 2Jfv /3 C/3 = , (288a) 
D( a - 2Mv' 3 C a ;/3 = 0. (288b) 
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They are equivalent to the homogeneous wave equations in the conformal flat spacetime 

□ «) = o , □ «") = o , 



(289) 



which point out that (in the approximation under consideration) the products, a( and a( a , are the harmonic functions. 

Potentials q,p,p a ,P a p must satisfy the gauge conditions (259a), (259b). Neglecting terms being quadratic with 
respect to the Hubble parameter, the gauge conditions (259a), (259b) can be written down as follows 



(ap a ) a +v a (2aq-a P ) a +Xap = 
(ap^) j j +v l3 (ap a )j3 + ^ a ^{ap)^+^ap a = 



(290a) 
(290b) 



where the potentials p a and p TQ/3 are obtained from p a and p J afj by rising the indices with the Minkowski metric and 
taking into account that the indices of r Q and rLg in the integrands of (285c) and (285d) should be raised with the 
full background metric g a ^ = a~ 2 rf^ taken at the point of integration. It yields 



ap a (t, x) 
ap Tal3 (t,x) 



■J 

Jv 

■J 

Jv 



a 4 [q{s,x')} T a (s,x')d 3 x' 

\x — x'\ 
a 5 [r)(s, x 1 )] (s,x')d 3 x' 



(291a) 
(291b) 



It is instrumental to write down solutions for the products of the potentials p and p a = 7] a ^pp with the Hubble 
parameter. Multiplying both sides of equations (282b), (282c) with the Hubble parameter Jf, and neglecting the 
quadratic with respect to 5f terms, we obtain 



U{%ap) = 167ra 3 JfT 
which solutions are the retarded potential 



□ (Jfop a ) = 16ira 4 j<T a 



jiap(t,x) = -4 / 
Jv 

- 4 I 



3<ap a (t,x) 



a 3 [77(5, x 1 )] 3i frjs, x')} t (s, x') d 3 x' 
v \x-x'\ 
' a 4 [rjjs, x')} Vi [r](s, x')} r (a, x') d 3 x' 
\x — x'\ 



(292) 

(293a) 
(293b) 



Substituting functions q,p,p a ,p Tal3 and 3iap, j~Cap a to the gauge equations (290a), (290b), bring about the following 
integral equations 



(a 4 T a +v a a 3 a) 



1 



a 3 3iT 



+ a 4 jir° 



\x — X' 



\x — X' 



0, 
0, 



(294a) 
(294b) 



where all functions in the integrands are taken at the retarded time s and at the point x' , for example, a — a[r)(s, x')], 
ji = ji[r](s,x')], g = a[(s,x')], and so on. These equations are satisfied by the equations of motion (105a), (105b) 
of the localized matter distribution. Indeed, divergences of any vector F a and a symmetric tensor F Q/3 obey the 
following equalities 



pa 



1 



(295) 
(296) 



Moreover, the root square of the determinant of the background metric tensor is expressed in terms of the scale 
factor, \f^-g = a 4 , while the four-velocity u a — v a /a. Applying these expressions along with equations (295), (296) 
in equations of motion (105a), (105b), transforms them to 



(a 4 T a + v a a 3 a) a + a 3 j<T 



(a 4 T^ + a 3 v^r a ) 



+ 2a 3 jiT° 



(297a) 
(297b) 
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Equation (297a) proves that the integral equation (294a) and, hence, the gauge condition (290a) are valid. In order 
to prove the second integral equation (294b), we multiply equation (297b) with the scale factor a, and reshuffle its 
terms. It brings (297b) to the following form 

(a 5 r a/3 + a 4 v T a ) + a A< K T a = . (298) 

Substituting, T ali = + (l/3a 2 )7f Q/3 T, to (298) and comparing with the integrand of (295) makes it clear that 
(294b) is valid. It proves the second gauge condition (290b). We conclude that the retarded integrals (285a)-(285d) 
yield the complete solution of the linearised wave equations (282a)-(282d). Thus, we can chose the gauge functions 

C = C" = o. 



3. The matter field perturbations 



What remains is to find out solutions for the scalar functions V m and V q and Xm and Xq- I n the linearized Hubble 
approximation equation for V m is obtained from (253) by discarding all terms of the order of !K 2 . It yields 



□Vm + ( X - ^) v Q v^V m ,^ + ^3 - J) MvT m ,„ = -Ana 2 (a + r) 



Applying relationships (250), (251) in equation (299) allows us to recast it to 



1 



□ (a n V m ) + ( 1 - % j v Q v^ (a n V m ) 



3+(2n-l)- 



-Ana 2 (a + r) 



where n is yet undetermined real number. Choosing, n — n s , with 



1 



n. = - 1 - 3-£ 



annihilates the term being proportional to JC in the left side of (300) and reduces it to 

□ K ls K0 + - J) v Q v^ (a n - V m ) >ap - -Ana 2+n ° (a + r) . 



(299) 



(300) 



(301) 



(302) 



This equation describes propagation of perturbation V m with the speed of sound c s . Indeed, let us introduce the 
sound-wave Laplace-Beltrami operator 



Then, equation (302) reads 



n s = ( --v a ^ + 7t a0 ) d a 



□s (fl" s V m ) - -47ra 2 +"= (a + r) 



(303) 



(304) 



This equation has a well-defined Green function with characteristics propagating with the speed of sound c s . We 
discard the advanced Green function because we assume that at infinity the function V m and its first derivatives 
vanish. Solution of (304) is explained in Appendix B, and has the following form 



a n °V m (t,x) = 



a 2+n °(^x')[<j(<;,x')+T(<;,x')] d 3 x' 



V Jl+(l-f) 7 2G9xn)S 



IX — X 



(305) 



where the retarded time <r is given by equation (B18), (3 — /3 l — v l /c, 7 = 1/yT — /3 2 is the Lorentz factor, and the 
unit vector n = (x — x')/\x — x'\. 

Linearized equation for Vq is obtained from (254) after discarding all terms being proportional to Jf 2 . It yields 



□V q + 2 (l - \[^^ = -4ira 2 (a + r) . 



(306) 
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Applying relationship (250) in (306) allows us to recast it to 



— □ (a n V q ) + 2\n + l-^ — Ax 2 j Kv a V nha = -Ana z (a + r) . (307) 

Choosing, n = n q = — 1 + ^/3/ (2xi)Ax 2 , eradicates the second term in the left side of (307) that results in 

□ (a"<%) = -47ra 2+rl " (a + r) . (308) 
This is the wave equation in flat spacetime. We pick up the retarded solution as the most physical one. It reads 

a 2+« q ( Sj x /) [(J + T ( S; jp/)] d 3 x , 



a n *V q 



I (s,x') [a( S ,xi) + T( S ,x')}d*x> ^ (30Q) 

Jv \ x ~ x '\ 

where the retarded time s has been defined in (286). 

Perturbations Xm and x q can be found by integrating equations (180a) and (180b) that can be written as 

v Q Xm,« = a (y m + |) , v Q x q ,a = a (V q + |) . (310) 
These are the ordinary differential equations of the first order. Their solutions are 

r* i 

X m = / a[t,x(t)]{V m [t,x(t)] + -q[t lX (t)]}dt, (311a) 

Jto 1 

X q - / a[t,x(t)]{V q [t,x{t)] + -q[t,x(t)]}dt , . (311b) 

where to is an initial epoch of integration, and the integration is performed along the characteristics of the unperturbed 
equations of motion of matter of the background universe 

^=v%x). (312) 

These characteristics make up the Hubble flow. Therefore, the most simple way to integrate equations (310) would be 
to work in the preferred coordinate frame X a = (ry,X 4 ) where the velocity v l — 0, and the coordinates X 1 = const. 
After the calculations in the rest frame of the Hubble flow are finished, the transformation to a moving frame can be 
done with the Lorentz boost. 
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Appendix A: Lorentz Invariance of the Retarded Potential 

We use a prime in the appendices exclusively as a notation for time and spatial coordinates which are used as 
variables of integration in volume integrals (see, for example, equations (A2), (A3), and so on). It should not be 
confused with the time derivative with respect to the conformal time used in the main text of the present paper. 

Let us consider an inhomogeneous wave equation for a scalar field, V = V(r/,X), written down in a coordinate 
chart X a = (X°, X 1 ) = (77, X), 

UV = -4iro- x , (Al) 

where □ = i] al3 d a p, d a — d/dX a , and a x = <J x (r\, X) is the source (a scalar function) of the field V with a compact 
support (bounded by a finite volume in space). Equation (Al) has a solution given as a linear combination of advanced 
and retarded potentials. Let us focus only on the retarded potential which is more common in physical applications. 
Advanced potential can be treated similarly. 

We assume the field, V, and its first derivatives vanish at past null infinity. Then, the retarded solution (retarded 
potential) of (Al) is given by an integral, 

L "% X X X ' - <A2) 

where 

( = r l -\X-X'\, (A3) 

is the retarded time, and we assume the fundamental speed c = 1. Physical meaning of the retardation is that the 
field V propagates in spacetime with the fundamental speed c from the source <r x , to the point with coordinates 
X a = (j],X) where the field V is measured in correspondence with equation (A2). Left side of equation (Al) is 
Lorentz- invariant. Hence, we expect that solution (A3) must be Lorentz-invariant as well. As a rule, textbooks prove 
this statement for a particular case of the retarded (Lienard-Wiechert) potential of a moving point-like source but 
not for the retarded potential given in the form of the integral (A2). This appendix fulfils this gap. 

Lorentz transformation to coordinates, x a — (t, x) linearly transforms the isotropic coordinates X a = (77, X) of the 
FLRW universe as follows 

/ = k a f)X 13 , (A4) 
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where the matrix of the Lorentz boost [79] 

A%=7, A* = A ^ = —7/3* , A*,- = ^ + ^/3*/3 j , (A5) 
the boost four-velocity u a — = u°{l,/3 4 } is constant, and 

7 = «° = -. 1 , (A6) 

is the constant Lorentz-factor. 

The inverse Lorentz transformation is given explicitly as follows 



77 = 7 (i + ■ x) , (A7) 
1 + 7 



X = r +---(/3T)6, (A8) 



where 

r = x + 0t , (A9) 

and the boost three-velocity, — {(3 1 } — {u l /u }. 

Let us reiterate (A2) by introducing a one-dimensional Dirac's delta function and integration with respect to time 

where Q is the retarded time given by (A3). Then, we transform coordinates X' a = (77', X') to x' a = (t 1 , x') with the 
Lorentz boost (A4). The Lorentz transformation changes functions entering the integrand of (A10) as follows, 

a(r)',X') = a x (t',x') , (All) 

\X-X'\ = v/|r-r'| 2 +7 2 [/3-(r-r')] 2 , (A12) 

where the coordinate difference 

r-r' = x-x' +0(t-t') . (A13) 

The coordinate volume of integration remains Lorentz-invariant 

dri'SX' = dt'cfx' . (A14) 

Let us denote F v (r]') = r]' — Q where ( is given by (A3). After making the Lorentz transformation this function changes 
to 



F„{rf) = F t (f) = 7 [f - t - ■ (x - x')} + y/\x - x'\ 2 - (t> - tf + 7 2 [/3 • (x - x>) - (t> - t)Y , (A15) 

where we have used equations (A7), (A8) and (A12) and relationship 7 2 /3 2 = 7 2 — 1, to perform the transformation. 
Integral (A10) in coordinates x a becomes 

v( t X )= T f °*v>*mw))<** i *' (A16) 

J -00 Jv y/\r - r'| 2 + 7 2 [/3 • (r - r')] 2 ' 

The delta function has a complicated argument F t {t') in coordinates x a . It can be simplified with a well-known 
formula 

6[F t (t')] = S X^-, (A17) 

where F t (s) = [dF t (t')/dt'] t , =s , and s is one of the roots of equation F t (t') — that is associated with the retarded 
interaction. It is straightforward to confirm by inspection that the root is given by formula, 

s = t-\x-x'\. (A18) 
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The time derivative of function F t (t') is 



HH)=lW , /3 2 (t' -t)-f3-(x-x') 

y/\x - x'\ 2 - (f - t) 2 + 1 2 [(3-{x- x>) - (F - t)] 2 

After substituting t' = s, with s taken from equation (A18), we obtain, 

" ^ |x-xl"l-x-) • < A20) 
Performing now integration with respect to t' in equation (A16) with the help of the delta-function, we arrive to 

v{t - x) = I F?*x' ] xn ■ < A21) 

Jv F t (s)\X - X'\ t , =s 

where \X — X'\ t > =s must be calculated from (A12) with t' = s where s is taken from (A18). It yields 

F t {s)\X - X'\ v=a = \x - x'\ , (A22) 
and proves that the retarded potential (A2) is Lorentz-invariant 

a x (C,X')d 3 A' f a x (s,x>) d s x> 



\X-X'\ J v \x-x'\ 

We have verified the Lorentz invariance for the scalar retarded potential. However, it is not difficult to check that 
it is valid in case of a source cr aia2 Q; that is a tensor field of rank /. Indeed, the Lorentz transformation of the 
source leads to A' 3l ctl A' 32 Q2 ...A^' Q; (T^ 1( 3 2 ... ( 3 1 but the matrix A a fj is constant, and can be taken out of the sign of the 
retarded integral. Because of this property, all mathematical operations given in this appendix for a scalar retarded 
potential, remain the same for the tensor of any rank. Hence, the Lorentz invariance of the retarded integral is a 
general property of the wave operator in Minkowski spacetime. 



Appendix B: Retarded Solution of the Sound- Wave Equation 

Let us consider an inhomogeneous sound-wave equation for a scalar field U = U(r], X) written down in the isotropic 
coordinates X a = (77, X), 

D S U = -Attt x , (Bl) 

where t x = T x (r], X) is the source of U having a compact support, and the sound-wave differential operator D s was 
defined in (303). It is Lorentz-invariant and reads 

□ B = □ + (l - V) V a v"0 a /9 , (B2) 



where v" is four- velocity of motion of the medium with respect to the coordinate chart, c s is the constant speed of 
sound in the medium, and we keep the fundamental speed c in the definition of the operator for dimensional purposes. 
We assume that c s < c. The case of c s = c is treated in section A, and the case of c s > c makes a formal mathematical 
sense in discussion of the speed of propagation of gravity in alternative theories of gravity since the equation describing 
propagation of gravitational potential U has the same structure as (Bl) after formal replacement of c s with the speed 
of gravity c g [59, 104]. In particular, in the Newtonian theory the speed of gravity c g = 00, and the operator (B2) is 
reduced to the Laplace operator 

A = □ + v a v /3 9 a a /3 = n af! d a p , (B3) 

where the constant projection operator, 7r a/3 , has been defined in (239). 

We are looking for the solution of (Bl) in the Cartesian coordinates x a = (t, x) moving with respect to the isotropic 
coordinates X a with constant velocity j3 l . Transformation from X a to x a is given by the Lorentz transformation 
(A4). In coordinates X a the four-velocity v Q = (1,0,0,0). Therefore, in these coordinates, equation (Bl) is just a 
wave equation for the field U propagating with speed c s . It has a well-known retarded solution, 
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where 

Vs = V --\X-X'\, (B5) 

Cs 

is the retarded time. 

Equation (Bl) is Lorentz- invariant. Hence, its solution must be Lorentz-invariant as well. Our goal is to prove this 
statement. To this end, we take solution (B4) and perform the Lorentz transformation (A7), (A8). We recast the 
retarded integral (B4) to another form with the help of one-dimensional delta-function 

It looks similar to (A2) but one has to remember that the retarded time rj s differs from ( that was defined in (A3) on 
the characteristics of the null cone defined by the fundamental speed c. Transformation of functions entering integrand 
in (B6) is similar to what we did in section A but, because c s ^ c, calculations become more involved. It turns out 
more preferable to handle the calculations in tensor notations, making transition to the coordinate language only at 
the end of the transformation procedure. 

Let us consider two events with the isotropic coordinates X a = (t],X) and X' a = (i]',X r ). We postulate that in 
the coordinate chart, x a , these two events have coordinates, x a = (t,x), and, x' a = (t',x'), respectively. We define 
the components of a four-vector, r a = (t' — t,x — x') which is convenient for doing mathematical manipulations with 
the Lorentz transformations. For instance, the Lorentz transformation of the Euclidean distance between the spatial 
coordinates of the two events, is given by a 

\X - X'\ = yJit a pr°rP , (B7) 

where 7f Q/3 is the operator of projection on the hyperplane being orthogonal to v Q (the same operator as in (B3)). 
Equation (B7) is a Lorentz-invariant analogue of expression (A12) and matches it exactly. Transformation of the 
source, T x {X a ) = r x (x a ) is fully equivalent to that of a x as given by equation (All). Coordinate volume of integration 
transforms in accordance with (A14). We need to transform the argument, rj' — r/s, of delta-function which we shall 
denote in coordinates X a as f v {r)') = rj' — r} s . The argument is a scalar function which is transformed as /^(r/) = ft(t') 
where, 

ft(t') = ~v a r a + -Jir aP r«rP . (B8) 
Transformation of the delta-function in the integrand of integral (B6) is 

5{Mt')} = 6 -^-^, (B9) 

/t(0 

where /*(<j) = [dft(t')/dt'] t , , and <; is one of the roots of equation /t(t') = that is associated with the retarded 
interaction. Eventually, after accounting for transformation of all functions and performing integration with respect 
to time, integral (B6) assumes the following form 

U(t,x)= j . ZkfrlfVjS! , (BIO) 

Jv Ms)\x - x>\ v= , 

where \X — X'| t / =s denotes the expression (B7) taken at the value of t' = <r. What remains is to calculate the instant 
of time, <r, and the value of functions entering denominator of the integrand in (BIO). 

Calculation of ? is performed by solving equation / t (<r) = 0, that defines the characteristic cone of the sound waves, 
and has the following explicit form, 



Va/3 + \ 1 - ^ ) v " v /3 



rV = , (Bll) 



which is derived from (B8). This is a quadratic algebraic equation with respect to the time variable r° = <r — t. It 
reads 

A{s-tf + 2B{q-t) + C = , (BI2) 
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where the coefficients A, B, C of the quadratic form are, 



A = -1+1- 



B = - 1 



7 2 /3 • (x - as') , 



C = \x-x'\ 2 + ( 1 



7 2 \J3 ■ (x - *')] . 



(B13) 
(B14) 
(B15) 



Equation (B12) has two roots corresponding to the advanced and retarded times. The root corresponding to the 
retarded-time solution of (B12) is 



or, more explicitly, 



? = t — \x — x 



? = t - 1 (B - VB 2 - AC) , 



1-1-37 



where the unit vector n = (a; — x')/\x — x'\. After some algebra equation (B17) can be simplified to 



a s 



c; = t \x — x , 



where 



1-/3 2 



1 - 



2 



'1+ 1 



7 2 (/3 x n) 2 



7 2 (/3-n) 



(B16) 



(B17) 



(B18) 



(B19) 



Coefficient a s defines the speed of propagation of the sound waves, v s = c s /a s , as measured by observer moving with 
speed ft 1 with respect to the Hubble flow. Thus, the value of the speed of sound, v B , depends crucially on the motion 
of observer. 

Derivative of the function, /*(<:), is given by 



where the partial derivative dr a /d<; = 6g = (1,0,0,0). Making use of (B8), the partial derivative 

dr a a c s ^fTt at) r a rP ' 



(B20) 



(B21) 



which has to be calculated at the instant of time, t' = <;, where q is given by (B18). 

In order to calculate the denominator in the integrand in (BIO), we account for (B7), (Bll) and combine (B20), 
(B21) together. We get 



\X - X'lU,) = - 



1-4 v Q v )9 r' : 



°a ■ 



(B22) 



It is straightforward to check that after using (B16) the above equation is reduced to \X— X'\f x {<;) = (c/c s )V B 2 — AC, 
or more explicitly, 



\X-X'\f x ^) = \x-x'\Jl+ 1 



7 2 (/3 x n) 2 



(B23) 
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Finally, the retarded Lorentz- invariant solution of (Bl) is 



r(Lx)= T X (,,X') dV (B24) 



'l+(l--|f(/3xn) 2 



with the retarded time <; calculated in accordance with (B18). This solution reduces to the retarded potential (A23) 
in the limit of c s — > c. 



